{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 变分量子奇异值分解\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览\n",
    "\n",
    "在本教程中，我们一起学习下经典奇异值分解（singular value decomposition, SVD）的概念以及我们自主研发的量子神经网络版本的量子奇异值分解（variational quantum SVD, VQSVD）[1] 是如何运作的。主体部分包括两个具体案例：\n",
    "- 分解随机生成的 $8\\times8$ 复数矩阵；\n",
    "- 应用在图像压缩上的效果。\n",
    "\n",
    "## 背景\n",
    "\n",
    "奇异值分解有非常多的应用，包括主成分分析（principal component analysis, PCA）、求解线性方程组和推荐系统。其主要任务是给定一个复数矩阵 $M \\in \\mathbb{C}^{m \\times n}$, 找到如下的分解形式：$M = UDV^\\dagger$。其中 $U_{m\\times m}$ 和 $V^\\dagger_{n\\times n}$ 是酉矩阵（Unitary matrix）, 满足性质 $UU^\\dagger = VV^\\dagger = I$。 \n",
    "\n",
    "- 矩阵 $U$ 的列向量 $|u_j\\rangle$ 被称为左奇异向量（left singular vectors）, $\\{|u_j\\rangle\\}_{j=1}^{m}$ 组成一组正交向量基。这些列向量本质上是矩阵 $MM^\\dagger$ 的本征向量。\n",
    "- 类似的，矩阵 $V$ 的列向量 $\\{|v_j\\rangle\\}_{j=1}^{n}$ 是 $M^\\dagger M$ 的本征向量也组成一组正交向量基。\n",
    "- 中间矩阵 $D_{m\\times n}$ 的对角元素上存储着由大到小排列的奇异值 $d_j$。 \n",
    "\n",
    "我们不妨先来看个简单的例子（为了方便讨论，我们假设以下出现的 $M$ 都是方阵）：\n",
    "\n",
    "$$\n",
    "M = 2*X\\otimes Z + 6*Z\\otimes X + 3*I\\otimes I = \n",
    "\\begin{bmatrix} \n",
    "3 &6 &2 &0 \\\\\n",
    "6 &3 &0 &-2 \\\\\n",
    "2 &0 &3 &-6 \\\\\n",
    "0 &-2 &-6 &3 \n",
    "\\end{bmatrix}, \\tag{1}\n",
    "$$\n",
    "\n",
    "那么该矩阵的奇异值分解可表示为:\n",
    "\n",
    "$$\n",
    "M = UDV^\\dagger = \n",
    "\\frac{1}{2}\n",
    "\\begin{bmatrix} \n",
    "-1 &-1 &1 &1 \\\\\n",
    "-1 &-1 &-1 &-1 \\\\\n",
    "-1 &1 &-1 &1 \\\\\n",
    "1 &-1 &-1 &1 \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix} \n",
    "11 &0 &0 &0 \\\\\n",
    "0 &7 &0 &0 \\\\\n",
    "0 &0 &5 &0 \\\\\n",
    "0 &0 &0 &1 \n",
    "\\end{bmatrix}\n",
    "\\frac{1}{2}\n",
    "\\begin{bmatrix} \n",
    "-1 &-1 &-1 &-1 \\\\\n",
    "-1 &-1 &1 &1 \\\\\n",
    "-1 &1 &1 &-1 \\\\\n",
    "1 &-1 &1 &-1 \n",
    "\\end{bmatrix}. \\tag{2}\n",
    "$$\n",
    "\n",
    "我们通过下面几行代码引入必要的 library和 package。\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.433321Z",
     "start_time": "2021-03-09T03:47:28.730090Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import pi as PI\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import unitary_group\n",
    "from scipy.linalg import norm\n",
    "\n",
    "import paddle\n",
    "from paddle_quantum.ansatz import Circuit\n",
    "from paddle_quantum.linalg import dagger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 画出优化过程中的学习曲线\n",
    "def loss_plot(loss):\n",
    "    '''\n",
    "    loss is a list, this function plots loss over iteration\n",
    "    '''\n",
    "    plt.plot(list(range(1, len(loss)+1)), loss)\n",
    "    plt.xlabel('iteration')\n",
    "    plt.ylabel('loss')\n",
    "    plt.title('Loss Over Iteration')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 经典奇异值分解\n",
    "\n",
    "那么在了解一些简单的数学背景之后， 我们来学习下如何用 Numpy 完成矩阵的奇异值分解。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.469286Z",
     "start_time": "2021-03-09T03:47:34.440399Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "我们想要分解的矩阵 M 是：\n",
      "[[ 3.+0.j  6.+0.j  2.+0.j  0.+0.j]\n",
      " [ 6.+0.j  3.+0.j  0.+0.j -2.+0.j]\n",
      " [ 2.+0.j  0.+0.j  3.+0.j -6.+0.j]\n",
      " [ 0.+0.j -2.+0.j -6.+0.j  3.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# 生成矩阵 M\n",
    "def M_generator():\n",
    "    I = np.array([[1, 0], [0, 1]])\n",
    "    Z = np.array([[1, 0], [0, -1]])\n",
    "    X = np.array([[0, 1], [1, 0]])\n",
    "    Y = np.array([[0, -1j], [1j, 0]])\n",
    "    M = 2 *np.kron(X, Z) + 6 * np.kron(Z, X) + 3 * np.kron(I, I)\n",
    "    return M.astype('complex64')\n",
    "\n",
    "print('我们想要分解的矩阵 M 是：')\n",
    "print(M_generator())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.541244Z",
     "start_time": "2021-03-09T03:47:34.489833Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "矩阵的奇异值从大到小分别是:\n",
      "[11.  7.  5.  1.]\n",
      "分解出的酉矩阵 U 是:\n",
      "[[-0.5+0.j -0.5+0.j  0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j -0.5+0.j -0.5+0.j -0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j -0.5+0.j  0.5+0.j]\n",
      " [ 0.5+0.j -0.5+0.j -0.5+0.j  0.5+0.j]]\n",
      "分解出的酉矩阵 V_dagger 是:\n",
      "[[-0.5+0.j -0.5+0.j -0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j -0.5+0.j  0.5+0.j -0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j  0.5+0.j  0.5+0.j]\n",
      " [-0.5+0.j  0.5+0.j -0.5+0.j -0.5+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# 我们只需要以下一行代码就可以完成 SVD \n",
    "U, D, V_dagger = np.linalg.svd(M_generator(), full_matrices=True)\n",
    "\n",
    "# 打印分解结果\n",
    "print(\"矩阵的奇异值从大到小分别是:\")\n",
    "print(D)\n",
    "print(\"分解出的酉矩阵 U 是:\")\n",
    "print(U)\n",
    "print(\"分解出的酉矩阵 V_dagger 是:\")\n",
    "print(V_dagger)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.600873Z",
     "start_time": "2021-03-09T03:47:34.565570Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.+0.j  6.+0.j  2.+0.j  0.+0.j]\n",
      " [ 6.+0.j  3.+0.j  0.+0.j -2.+0.j]\n",
      " [ 2.+0.j  0.+0.j  3.+0.j -6.+0.j]\n",
      " [ 0.+0.j -2.+0.j -6.+0.j  3.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# 再组装回去，能不能复原矩阵？\n",
    "M_reconst = np.matmul(U, np.matmul(np.diag(D), V_dagger))\n",
    "print(M_reconst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "那当然是可以复原成原来的矩阵 $M$ 的！读者也可以自行修改矩阵，试试看不是方阵的情况。\n",
    "\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 量子奇异值分解\n",
    "\n",
    "接下来我们来看看量子版本的奇异值分解是怎么一回事。简单的说，我们把矩阵分解这一问题巧妙的转换成了优化问题。通过以下四个步骤：\n",
    "\n",
    "- 准备一组正交向量基 $\\{|\\psi_j\\rangle\\}$, 不妨直接取计算基 $\\{ |000\\rangle, |001\\rangle,\\cdots |111\\rangle\\}$ （这是3量子比特的情形）\n",
    "- 准备两个参数化的量子神经网络 $U(\\theta)$ 和 $V(\\phi)$ 分别用来学习左/右奇异向量\n",
    "- 利用量子神经网络估算奇异值 $m_j = \\text{Re}\\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle$\n",
    "- 设计损失函数并且利用飞桨来优化\n",
    "\n",
    "$$\n",
    "L(\\theta,\\phi) = \\sum_{j=1}^T q_j\\times \\text{Re} \\langle\\psi_j|U(\\theta)^{\\dagger} M V(\\phi)|\\psi_j\\rangle, \\tag{3}\n",
    "$$\n",
    "\n",
    "其中 $q_1>\\cdots>q_T>0$ 是可以调节的权重（超参数）, $T$ 表示我们想要学习到的阶数（rank）或者可以解释为总共要学习得到的奇异值个数。\n",
    "\n",
    "\n",
    "\n",
    "### 案例1：分解随机生成的 $8\\times8$ 复数矩阵\n",
    "\n",
    "接着我们来看一个具体的例子，这可以更好的解释整体流程。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.656799Z",
     "start_time": "2021-03-09T03:47:34.620380Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "我们想要分解的矩阵 M 是：\n",
      "[[6.+1.j 3.+9.j 7.+3.j 4.+7.j 6.+6.j 9.+8.j 2.+7.j 6.+4.j]\n",
      " [7.+1.j 4.+4.j 3.+7.j 7.+9.j 7.+8.j 2.+8.j 5.+0.j 4.+8.j]\n",
      " [1.+6.j 7.+8.j 5.+7.j 1.+0.j 4.+7.j 0.+7.j 9.+2.j 5.+0.j]\n",
      " [8.+7.j 0.+2.j 9.+2.j 2.+0.j 6.+4.j 3.+9.j 8.+6.j 2.+9.j]\n",
      " [4.+8.j 2.+6.j 6.+8.j 4.+7.j 8.+1.j 6.+0.j 1.+6.j 3.+6.j]\n",
      " [8.+7.j 1.+4.j 9.+2.j 8.+7.j 9.+5.j 4.+2.j 1.+0.j 3.+2.j]\n",
      " [6.+4.j 7.+2.j 2.+0.j 0.+4.j 3.+9.j 1.+6.j 7.+6.j 3.+8.j]\n",
      " [1.+9.j 5.+9.j 5.+2.j 9.+6.j 3.+0.j 5.+3.j 1.+3.j 9.+4.j]]\n",
      "矩阵的奇异值从大到小分别是:\n",
      "[54.83484985 19.18141073 14.98866247 11.61419557 10.15927045  7.60223249\n",
      "  5.81040539  3.30116001]\n"
     ]
    }
   ],
   "source": [
    "# 先固定随机种子， 为了能够复现结果\n",
    "np.random.seed(42)\n",
    "\n",
    "# 设置量子比特数量，确定希尔伯特空间的维度\n",
    "N = 3\n",
    "\n",
    "# 制作随机矩阵生成器\n",
    "def random_M_generator():\n",
    "    M = np.random.randint(10, size = (2**N, 2**N)) + 1j*np.random.randint(10, size = (2**N, 2**N))\n",
    "    M1 = np.random.randint(10, size = (2**N, 2**N)) \n",
    "    return M\n",
    "\n",
    "M = random_M_generator()\n",
    "M_err = np.copy(M)\n",
    "\n",
    "# 打印结果\n",
    "print('我们想要分解的矩阵 M 是：')\n",
    "print(M)\n",
    "\n",
    "U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
    "print(\"矩阵的奇异值从大到小分别是:\")\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.692350Z",
     "start_time": "2021-03-09T03:47:34.671114Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "选取的等差权重为：\n",
      "[24.+0.j 21.+0.j 18.+0.j 15.+0.j 12.+0.j  9.+0.j  6.+0.j  3.+0.j]\n"
     ]
    }
   ],
   "source": [
    "# 超参数设置\n",
    "N = 3           # 量子比特数量\n",
    "T = 8           # 设置想要学习的阶数\n",
    "ITR = 100       # 迭代次数\n",
    "LR = 0.02       # 学习速率\n",
    "SEED = 14       # 随机数种子\n",
    "\n",
    "# 设置等差的学习权重\n",
    "weight = np.arange(3 * T, 0, -3).astype('complex128')\n",
    "print('选取的等差权重为：')\n",
    "print(weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们搭建如下的量子神经网络结构："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:47:34.725007Z",
     "start_time": "2021-03-09T03:47:34.702861Z"
    }
   },
   "outputs": [],
   "source": [
    "# 设置电路参数\n",
    "cir_depth = 20                           # 电路深度\n",
    "block_len = 2                            # 每个模组的长度"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义量子神经网络\n",
    "def U_theta() -> Circuit:\n",
    "\n",
    "    # 用 UAnsatz 初始化网络\n",
    "    cir = Circuit(N)\n",
    "    \n",
    "    # 搭建层级结构：\n",
    "    for _ in range(cir_depth):\n",
    "        cir.ry()\n",
    "        cir.rz()\n",
    "        cir.cnot()\n",
    "\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们来完成算法的主体部分：\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:49:24.011232Z",
     "start_time": "2021-03-09T03:47:40.712257Z"
    }
   },
   "outputs": [],
   "source": [
    "class NET(paddle.nn.Layer):\n",
    "    def __init__(self, matrix: np.ndarray, weights: np.ndarray) -> None:\n",
    "        super(NET, self).__init__()\n",
    "        \n",
    "        # 创建用来学习 U_dagger 的量子神经网络\n",
    "        self.cir_U = U_theta()\n",
    "        \n",
    "        # 创建用来学习 V_dagger 的量子神经网络\n",
    "        self.cir_V = U_theta()\n",
    "        \n",
    "        # 将 Numpy array 转换成 Paddle 支持的 Tensor\n",
    "        self.M = paddle.to_tensor(matrix)\n",
    "        self.weight = paddle.to_tensor(weights)\n",
    "\n",
    "    # 定义损失函数和前向传播机制\n",
    "    def forward(self):\n",
    "        \n",
    "        # 获取量子神经网络的酉矩阵表示\n",
    "        U = self.cir_U.unitary_matrix()\n",
    "        U_dagger = dagger(U)\n",
    "        \n",
    "        V = self.cir_V.unitary_matrix()\n",
    "        V_dagger = dagger(V)\n",
    "        \n",
    "        # 初始化损失函数和奇异值存储器\n",
    "        loss = 0 \n",
    "        singular_values = np.zeros(T)\n",
    "        \n",
    "        # 定义损失函数\n",
    "        for i in range(T):\n",
    "            loss -= paddle.real(self.weight)[i] * paddle.real(U_dagger @ self.M @ V)[i][i]\n",
    "            singular_values[i] = paddle.real(U_dagger @ self.M @ V)[i][i].numpy()\n",
    "        \n",
    "        # 函数返回两个矩阵 U 和 V_dagger 学习的奇异值以及损失函数    \n",
    "        return U, V_dagger, loss, singular_values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: -262.8009\n",
      "iter: 10 loss: -1886.9825\n",
      "iter: 20 loss: -2182.6929\n",
      "iter: 30 loss: -2280.7289\n",
      "iter: 40 loss: -2315.2391\n",
      "iter: 50 loss: -2333.0160\n",
      "iter: 60 loss: -2342.7280\n",
      "iter: 70 loss: -2348.7263\n",
      "iter: 80 loss: -2353.2796\n",
      "iter: 90 loss: -2356.9237\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkk0lEQVR4nO3deZhcdZ3v8feneu+k053QnZANkkgAA5fNEIEBx4so6KC4oKCj4jIyemU2cRTUe9V7h/ugjlfHUVEUFNwQQRSRAUURQdkSIIGAkRASkpClE9LZ0+v3/nFOJ5WmE7u7uvp0V31ez1NPVf3OqTrf0yfpT5/f7yyKCMzMzAqRy7oAMzMb+xwmZmZWMIeJmZkVzGFiZmYFc5iYmVnBHCZmZlYwh4mZDYqkHZLmZF2HjS4OExtzJK2UdFZGyz5N0m8lbZe0VdIvJM0bweXvXXdJ75F0X5GX9ztJf5ffFhHjI2JFMZdrY4/DxGyAJJ0K/Ar4OTANmA0sBv4w3H+pK1HU/5+SKov5/VZeHCZWMiTVSPqypOfTx5cl1aTTmiXdJqlN0guS7u39ZS3p45LWpnsbyyS96gCL+DxwfUT8R0Rsj4gXIuJTwAPAZ9LvekrSuXk1VUpqlXRS+v4USX9M61gs6ZV58/5O0hWS/gDsAg4YUJJeCnwDODXtdmrL+xn8u6TnJG2Q9A1Jdem0V0pak67veuA7kiamP5dWSVvS1zPS+a8AzgC+mi7jq2l7SDoifd0o6fr086skfSrv5/oeSfel9WyR9Kyk1w54g9qY4jCxUvJJ4BTgBOB4YAHwqXTapcAaoAWYAnwCCElHAZcAJ0dEA3A2sLLvF0uqB04DftLPcm8EXp2+/hHw9rxpZwObIuIRSdOBXwL/BkwCPgrcLKklb/53ARcDDcCqA61oRDwFfBC4P+12akonXQkcmf4MjgCmA/8r76OHpss+PF1ODvhO+v4wYDfw1XQZnwTuBS5Jl3FJP6X8J9BIEnx/DbwbeG/e9JcDy4BmkjC+RpIOtF42djlMrJT8LfC/I2JjRLQCnyX55QzQCUwFDo+Izoi4N5IL03UDNcA8SVURsTIinunnuyeR/H9Z18+0dSS/LAF+CLwhDR+Ad5AEDMA7gdsj4vaI6ImIXwMLgdflfdd3I2JpRHRFROdgVj79JX0x8C/pXtN24P8CF+bN1gN8OiLaI2J3RGyOiJsjYlc6/xUkoTCQ5VWk3315uqe2Evgi+37mAKsi4lsR0Q1cR7INpgxmvWxscJhYKZnG/n/Nr0rbAL4ALAd+JWmFpMsAImI58M8k3VQbJd0gaRovtoXkF/HUfqZNBTblfd9TwOvTQHkDScBA8tf/W9Murra0a+r0Pt+5ejAr3EcLUA8syvv+O9L2Xq0Rsaf3jaR6Sd9Mu6i2Ab8HmtKg+EuagSpe/DOfnvd+fe+LiNiVvhw/iHWyMcJhYqXkeZJf2L0OS9tI/3K+NCLmkPyC/0jv2EhE/DAiTk8/G8Dn+n5xROwE7gfe2s9y3wb8Ju99b1fXecCTacBAEhTfi4imvMe4iLgyf1GDWN++824i6aY6Ju/7GyNi/EE+cylwFPDyiJgAvCJt1wHm77u8Tl78M187iHWwEuEwsbGqSlJt3qOS5Jf4pyS1SGomGSv4PoCkcyUdkXYFbSXp3uqRdJSkM9OB+j0kv4x7DrDMy4CLJP2jpIZ08PrfgFNJutR63QC8BvgQ+/ZKSGt5vaSzJVWkdb+yd8B7CDYAMyRVA0RED/At4EuSJqfrPV3S2Qf5jgaSdW6TNAn4dD/L6PdAgLTr6kbgivTncTjwkXQ9rcw4TGysup3kl2Dv4zMkA9sLgSXA48AjaRvAXOAuYAfJHsbXI+JukvGSK0n+yl4PTAYu72+BEXEfyYD6m0nGSVYBJwKnR8TTefOtS5dxGvDjvPbVJHsrnwBaSfZU/pWh/z/8LbAUWC9pU9r2cZLuvAfSbqu7SPY8DuTLQB3J+j9A0i2W7z+A89Ojsb7Sz+f/AdgJrADuIwnPa4e0NjamyTfHMjOzQnnPxMzMCuYwMTOzgjlMzMysYA4TMzMrWNle6K25uTlmzZqVdRlmZmPKokWLNkVES9/2sg2TWbNmsXDhwqzLMDMbUyT1e804d3OZmVnBHCZmZlYwh4mZmRXMYWJmZgVzmJiZWcEcJmZmVjCHiZmZFcxhMki3PLqGHzx4wFtzm5mVJYfJIP1yyTp+8MBzWZdhZjaqOEwGqbGumq27O7Muw8xsVHGYDFJTfRVtuzqyLsPMbFRxmAxSU10VOzu66eg60G3CzczKj8NkkJrqqwDc1WVmlsdhMkiN9dWAw8TMLJ/DZJAa63r3TDxuYmbWy2EySE1pmLTt8p6JmVkvh8kg9Y6ZOEzMzPZxmAxSU10yZtLmMRMzs70cJoPUUFuJBFt9romZ2V4Ok0HK5cSE2iofzWVmlsdhMgRN9VXu5jIzy+MwGYKmuioPwJuZ5XGYDEFjfbX3TMzM8jhMhqCprsoD8GZmeRwmQ9BU7wF4M7N8DpMhaKxLwqSnJ7IuxcxsVHCYDEFjXRU9Advbu7IuxcxsVBh1YSLpM5LWSnosfbwub9rlkpZLWibp7Lz2c9K25ZIuK3aNTb1XDvYRXWZmAFRmXcABfCki/j2/QdI84ELgGGAacJekI9PJXwNeDawBHpZ0a0Q8Wazi9l7scXcHh1FfrMWYmY0ZozVM+nMecENEtAPPSloOLEinLY+IFQCSbkjnLV6Y+AZZZmb7GXXdXKlLJC2RdK2kiWnbdGB13jxr0rYDtReNrxxsZra/TMJE0l2SnujncR5wFfAS4ARgHfDFYVzuxZIWSlrY2to65O+ZsLeby2FiZgYZdXNFxFkDmU/St4Db0rdrgZl5k2ekbRykve9yrwauBpg/f/6Qj+vde7dFn7hoZgaMwm4uSVPz3r4JeCJ9fStwoaQaSbOBucBDwMPAXEmzJVWTDNLfWswaayorqK+ucDeXmVlqNA7Af17SCUAAK4G/B4iIpZJuJBlY7wI+HBHdAJIuAe4EKoBrI2JpsYtsqvOVg83Meo26MImIdx1k2hXAFf203w7cXsy6+mqsr/bRXGZmqVHXzTVWJBd7dJiYmYHDZMga66po2+0BeDMzcJgMWVO9b5BlZtbLYTJEjb51r5nZXg6TIWqqq6ajq4c9nd1Zl2JmljmHyRD5kipmZvs4TIaoMe/KwWZm5c5hMkR7L0PvPRMzM4fJUDW6m8vMbC+HyRD13m1xm4/oMjNzmAxVk8dMzMz2cpgMUX11BVUVcjeXmRkOkyGTlF5SxWFiZuYwKUCjL/ZoZgY4TArSVF/tMRMzMxwmBWmqq/I9TczMcJgUpKm+mi07HSZmZg6TArQ01NC6vZ2IyLoUM7NMOUwKMLmhho7uHh8ebGZlz2FSgMkTagDYuL0940rMzLLlMCnA5IZaADZu35NxJWZm2XKYFGByQ7pnss17JmZW3hwmBXA3l5lZwmFSgPrqSsbXVLqby8zKnsOkQJMbarxnYmZlz2FSoJaGGlo9ZmJmZc5hUqDJE2rdzWVmZc9hUiB3c5mZOUwKNrmhhl0d3exo78q6FDOzzDhMCtR7ePCGbe7qMrPy5TAp0N6z4D0Ib2ZlzGFSoL1nwXsQ3szKmMOkQL17Jq0ehDezMuYwKdCEukqqK3M+osvMyprDpECSmDKhho0egDezMuYwGQaTG2q9Z2JmZc1hMgx84qKZlTuHyTCY3OBuLjMrbw6TYTB5Qi3b9nSxp7M761LMzDKRSZhIequkpZJ6JM3vM+1yScslLZN0dl77OWnbckmX5bXPlvRg2v5jSdUjuS6QXDkYfHiwmZWvrPZMngDeDPw+v1HSPOBC4BjgHODrkiokVQBfA14LzAPens4L8DngSxFxBLAFeP/IrMI+PnHRzMpdJmESEU9FxLJ+Jp0H3BAR7RHxLLAcWJA+lkfEiojoAG4AzpMk4EzgpvTz1wFvLPoK9OFLqphZuRttYybTgdV579ekbQdqPwRoi4iuPu39knSxpIWSFra2tg5b0b7Yo5mVu8pifbGku4BD+5n0yYj4ebGWezARcTVwNcD8+fNjuL53Un01lTn58GAzK1tFC5OIOGsIH1sLzMx7PyNt4wDtm4EmSZXp3kn+/CMmlxPN432uiZmVr9HWzXUrcKGkGkmzgbnAQ8DDwNz0yK1qkkH6WyMigLuB89PPXwRkstczeYLDxMzKV1aHBr9J0hrgVOCXku4EiIilwI3Ak8AdwIcjojvd67gEuBN4CrgxnRfg48BHJC0nGUO5ZmTXJuETF82snBWtm+tgIuIW4JYDTLsCuKKf9tuB2/tpX0FytFemWhpqeeS5tqzLMDPLxGjr5hqzZk6q44WdHWzf05l1KWZmI85hMkzmNI8H4NlNOzOuxMxs5DlMhslLWsYBsKLVYWJm5cdhMkwOO6SenGBF646sSzEzG3EOk2FSU1nBzEn1PONuLjMrQw6TYTSneZy7ucysLDlMhtGclvE8u2kHPT3DdqUWM7MxwWEyjOa0jGNPZw/rfPKimZUZh8kw6j082IPwZlZuHCbDyIcHm1m5cpgMo5aGGsbXVHrPxMzKjsNkGEliTss4VvjwYDMrMw6TYebDg82sHDlMhtmclvGsbdvN7o7urEsxMxsxDpNhNicdhPcFH82snDhMhtnew4M3eRDezMqHw2SYzW724cFmVn4GFCaS/knSBCWukfSIpNcUu7ixqK66gulNdT482MzKykD3TN4XEduA1wATgXcBVxatqjFudrMPDzaz8jLQMFH6/DrgexGxNK/N+pjTkhweHOELPppZeRhomCyS9CuSMLlTUgPQU7yyxrajDm1gR3sXa7bszroUM7MRUTnA+d4PnACsiIhdkiYB7y1aVWPc8TOaAHhsdRszJ9VnW4yZ2QgY6J7JqcCyiGiT9E7gU8DW4pU1th11aAPVlTmWrGnLuhQzsxEx0DC5Ctgl6XjgUuAZ4PqiVTXGVVXkmDd1AovXOG/NrDwMNEy6IhlNPg/4akR8DWgoXllj3/EzGnli7Va6fddFMysDAw2T7ZIuJzkk+JeSckBV8coa+46b0cSujm6e8fkmZlYGBhomFwDtJOebrAdmAF8oWlUl4PiZTUAyCG9mVuoGFCZpgPwAaJR0LrAnIjxmchBzmsfRUFPpQXgzKwsDvZzK24CHgLcCbwMelHR+MQsb63I5cez0RpZ4EN7MysBAzzP5JHByRGwEkNQC3AXcVKzCSsFxMxu59r5nae/qpqayIutyzMyKZqBjJrneIEltHsRny9bxM5ro7A7+tG571qWYmRXVQPdM7pB0J/Cj9P0FwO3FKal09A7CL17Ttve1mVkpGlCYRMS/SnoL8Fdp09URcUvxyioN0xpraR5fzeLVW5NrCJiZlaiB7pkQETcDNxexlpIjieNmNPmILjMreQcNE0nbgf5O4RYQETGhKFWVkONmNHL3so1s39NJQ63P8zSz0nTQQfSIaIiICf08GhwkA/OywycSAY8815Z1KWZmReMjsorspMMmUpETDz/7QtalmJkVjcOkyMbVVHLs9EYecpiYWQnLJEwkvVXSUkk9kubntc+StFvSY+njG3nTXibpcUnLJX1FktL2SZJ+Lenp9HliFut0MC+fPYnHVrexp7M761LMzIoiqz2TJ4A3A7/vZ9ozEXFC+vhgXvtVwAeAuenjnLT9MuA3ETEX+E36flRZMGsSHd09LPZFH82sRGUSJhHxVEQsG+j8kqYCEyLigfS+KtcDb0wnnwdcl76+Lq991Dh51iQk3NVlZiVrNI6ZzJb0qKR7JJ2Rtk0H1uTNsyZtA5gSEevS1+uBKQf6YkkXS1ooaWFra+uwF34gjfVVHDWlgYdWOkzMrDQVLUwk3SXpiX4e5x3kY+uAwyLiROAjwA8lDfgQ5HSv5YC3NoyIqyNifkTMb2lpGfC6DIcFsyexaNUWOrt7RnS5ZmYjYcBnwA9WRJw1hM+0k9yEi4hYJOkZ4EhgLckNuXrNSNsANkiaGhHr0u6w/AtSjhoLZk/i+vtXsfT5bZzg63SZWYkZVd1cklokVaSv55AMtK9Iu7G2STolPYrr3cDP04/dClyUvr4or31UWTBrEgAPPbs540rMzIZfVocGv0nSGpLLH/4yvSIxwCuAJZIeI7lXygcjoneg4X8A3waWA88A/5W2Xwm8WtLTwFnp+1Fn8oRaZjeP46Fnt2RdipnZsCtaN9fBpFccftFVhw92McmIWAgc20/7ZuBVw11jMSyYNYk7lq6npyfI5ZR1OWZmw2ZUdXOVugWzJ7F1dydPrtuWdSlmZsPKYTKCzjiyGYDfLRuVxwiYmQ2Zw2QETW6o5bgZjfz2Tw4TMystDpMRdubRk3l0dRubd7RnXYqZ2bBxmIywM4+eTATc8+eROwPfzKzYHCYj7NhpjTSPr3FXl5mVFIfJCMvlxJlHt3DPn1t9aRUzKxkOkwycefRktu/pYtEqn8BoZqXBYZKB0+e2UFUh7nZXl5mVCIdJBsbXVPLy2YfwG4eJmZUIh0lG/vvRk1m+cQerNu/MuhQzs4I5TDLy2mMPJSe44eHVWZdiZlYwh0lGpjXV8ep5U7jhoefY09mddTlmZgVxmGTootNmsWVXJ7cufj7rUszMCuIwydCpcw7hqCkNXPfHlSR3HDYzG5scJhmSxLtPO5ylz2/zOSdmNqY5TDL2xhOm01BbyXf/uDLrUszMhsxhkrFxNZW8bf5M7nhiPeu37sm6HDOzIXGYjAIXnTqLnMRnbl3qsRMzG5McJqPAYYfU85HXHMkdS9fzs8fWZl2OmdmgOUxGiQ+cMYf5h0/kf/18Kc+37c66HDOzQXGYjBIVOfHFtx1Pd0/wsZuW0NPj7i4zGzscJqPI4YeM45N/81LuW76Jq+55JutyzMwGzGEyyrxjwWG8/vhpfOHOZfzCZ8ab2RhRmXUBtj9JfOH841jXtptLf7KYqY21zJ81KeuyzMwOynsmo1BtVQVXv3s+0xpr+cD1C32ZejMb9Rwmo9SkcdV8570LCOBD33/EVxY2s1HNYTKKzW4exxffejxPrtvGv/3yyazLMTM7IIfJKPeql07h718xh+8/8By3LfGAvJmNTg6TMeCjZx/FSYc1cdnNj7Nyk8dPzGz0cZiMAVUVOf7zHSdRkROX/mSxT2g0s1HHYTJGTG+q49Ovn8eiVVv43gOrsi7HzGw/DpMx5E0nTuevj2zhc3f8iTVbdmVdjpnZXg6TMUQSV7zpWARc/tPHfbl6Mxs1HCZjzIyJ9Xz8tUdz79ObuGnRmqzLMTMDHCZj0jtffjgnz5rI/7ntSTZs890ZzSx7DpMxKJcTnz//eDq6e9zdZWajgsNkjJrdPI5/Pftofvunjfz0Ed+d0cyylUmYSPqCpD9JWiLpFklNedMul7Rc0jJJZ+e1n5O2LZd0WV77bEkPpu0/llQ9wquTmfeeNouTZ03ks79Y6u4uM8tUVnsmvwaOjYjjgD8DlwNImgdcCBwDnAN8XVKFpArga8BrgXnA29N5AT4HfCkijgC2AO8f0TXJUH5310d/sphun8xoZhnJJEwi4lcR0ZW+fQCYkb4+D7ghItoj4llgObAgfSyPiBUR0QHcAJwnScCZwE3p568D3jhCqzEqzG4ex6dffwz3Pr2JL/36z1mXY2ZlajSMmbwP+K/09XRgdd60NWnbgdoPAdrygqm3vV+SLpa0UNLC1tbWYSo/e29fcBgXzJ/JV+9ezq+Wrs+6HDMrQ0ULE0l3SXqin8d5efN8EugCflCsOvJFxNURMT8i5re0tIzEIkfMZ887huNmNHLpjYtZ0boj63LMrMwULUwi4qyIOLafx88BJL0HOBf429h3bOtaYGbe18xI2w7UvhloklTZp73s1FZVcNU7X0ZVZY73X7eQ1u3tWZdkZmUkq6O5zgE+BrwhIvIvMnUrcKGkGkmzgbnAQ8DDwNz0yK1qkkH6W9MQuhs4P/38RcDPR2o9RpvpTXV8690vY/3WPbzrmgdp29WRdUlmViayGjP5KtAA/FrSY5K+ARARS4EbgSeBO4APR0R3OiZyCXAn8BRwYzovwMeBj0haTjKGcs3Irsro8rLDJ/Gtd89nRetOLrr2Ibbv6cy6JDMrAyrXs6fnz58fCxcuzLqMornryQ188PuLOGFmE9e852Qa66qyLsnMSoCkRRExv2/7aDiay4rgrHlT+MrbT2TxmjYu+Ob9PqnRzIrKYVLCXvffpvKd9yxg9Qu7eMtVf/RRXmZWNA6TEnf63GZ+dPEp7O7o5i1X/ZF7ny6d82vMbPRwmJSB42Y0cfOHTmNyQy0XXfsQX//dcl9p2MyGlcOkTMxqHsctHz6NvzluGp+/YxkXf28RG7d7HMXMhofDpIzUV1fylQtP4H+eO497lrVy1hfv4YcPPkePLxBpZgVymJQZSbz/9Nnc8c9nMG/aBD5xy+O89Zv3c/8zm7MuzczGMIdJmZrTMp4ffeAUPn/+cax+YRdv/9YDvO2b93Pf05s8nmJmg+aTFo09nd3c8NBzXHXPM2zY1s68qRN43+mzef3xU6mprMi6PDMbRQ500qLDxPba09nNzx5dy7V/eJY/b9hB8/hq3nnK4bzzlMNpHl+TdXlmNgo4TPpwmBxYRPCH5Zv59n0r+N2yVqorc7zphOm8/4zZHDmlIevyzCxDBwqTyv5mtvImidPnNnP63GaWb9zOtX9YyU8fWcOPF67mr49s4e/OmM3pRzST3OjSzMx7JlmXMWa8sLODHz64iuvuX0Xr9nZe0jKOd51yOG9+2Qwm1Poikmblwt1cfThMhqa9q5tfLlnH9fev4rHVbdRXV/D646ZxwYKZnDizyXsrZiXOYdKHw6Rwj6/ZyvceWMkvFq9jd2c3R04ZzwUnH8abTpzOpHHVWZdnZkXgMOnDYTJ8tu/p5LYl67jh4dUsXt1GdUWOV8+bwgUnz+T0I5rJ5by3YlYqHCZ9OEyK40/rt/Hjh1dzy6NradvVyfSmOs5/2QzectIMDjukPuvyzKxADpM+HCbF1d7Vza+f3MCPH17Nfcs3EQHHz2jk3OOm8ZpjpnDYpHqPr5iNQQ6TPhwmI2fNll3ctmQdty15nifWbgPg0Am1vHzOJE6c2cSclvG8ZPJ4pk6odZeY2SjnMOnDYZKNlZt2cu/yTTy4YjMPrHiBTTva906rrshxaGMtUxtrmd5Ux4yJdcyYWM+MiXUc3jzOYWM2CvikRRsVZjWPY1Zzco5KRNC6o51nNu5kxaYdPLd5F89v3cO6tt3cv2IzG7btIf/q+DWVOQ4/pJ7DJtUzc1I9MyfWc2hjLZMbapjcUEvTuCoaairdfWaWAYeJZUYSkxtqmdxQy6kvOeRF0zu7e1i/dQ/PvbCLlZt3snLTTlZu3sXqF3bxx2c2s6uj+0WfqciJCbWVNNZVMaGuigm1VUwcV82k+iomjath8oQaDm2s5dAJtUxrqqOxzidcmg0Hh4mNWlUVuWQPZFI9f3VE837TIoIXdnawcXs7G7btoXV7O1t3d9K2q5O23R1s293Ftj2dbN3dydq23Wze0c62PV0vWkZDTSXT87rTZkysY3pTHdOa6pjaVEvzuBp3rZkNgMPExiRJHDK+hkPG1/DSqRMG9JmOrh5ad7Szfutu1m3dw/Ntu1m7ZTert+xm9Qu7uP+ZTezss7dTVZHsPU2ZkHSlTZ5QQ8v4GpobajhkXPXe56b6ahpqKh08VrYcJlY2qitzTG9K9jz6ExG07Ur2ZNZt3cO6NHQ2bEsey1t3cP+KzWzd3dnv53OCxroqxtVUMr6mknE1ldRW5aitrKCmKkdVxb5HdYWS58rc3ueayr7PFdSkz7VV+55rqyrSR466qgoqK3yPO8uew8QsJYmJ46qZOK6aY6c3HnC+9q5uNu3o4IUdHWza2c7mHR207eqgbVfSrbazvYsd7V3s7OiivbOHrbs7ae/sobO7h87uoL2r93UPHV09dPUUdkRlRU5p6CSBU52GUXUaUnsDqmJfcOXPU1O5f3t+2FVV5KisEJW55Lmq93VOVKbTqvpOSz9bkUumVVQomT8nKnLyARIlymFiNkg1lRUH3cMZrJ6eoLMnCZaOrh7au3rY09lNR3cP7Z3J6962Pelze2c3e140rXu/7+jsTp47unrY0d61d1pviHXkTe/o7mGkzhLoDZXe597gqcyJijScKnrnqRAVaXjlf6b3dU7J61zvNO2b3vvIaf/P5bfletuUvhZ7v29fm/Zrq8ix33Jz6p2XF82b623P+77e+ZUuqyJ/PcTeeXIifd43bTQHscPELGO5nKjJVWR+i+SudM+po6uHzp7995y6uoPO7uR1d08yX1d3EoJd3XltedO60mldPfu+o/fz+14n39sTsV9bb3t3T9LW1dOzt213Z9DTs/+83bHvdRLOyXNv+37zR4xYcA43if1CqjdsDhRCve29gaV02rUXnTzslzdymJgZQNptBXXV2YbaSNgvaCIJmuhhv7a94bRfG/tPjyAG0d4T+76nJ/YFXk9P0BP75o1Iatk7T+/8eZ+PPt8VkX5XkC4773t6goC9bTVVwz/O5jAxs7KTy4kcoqr0c3PE+DAQMzMrmMPEzMwK5jAxM7OCOUzMzKxgDhMzMyuYw8TMzArmMDEzs4I5TMzMrGBle9teSa3AqkF8pBnYVKRyRqtyXGcoz/Uux3WG8lzvQtf58Iho6dtYtmEyWJIW9nff41JWjusM5bne5bjOUJ7rXax1djeXmZkVzGFiZmYFc5gM3NVZF5CBclxnKM/1Lsd1hvJc76Kss8dMzMysYN4zMTOzgjlMzMysYA6Tv0DSOZKWSVou6bKs6ykWSTMl3S3pSUlLJf1T2j5J0q8lPZ0+T8y61uEmqULSo5JuS9/PlvRgus1/LKk66xqHm6QmSTdJ+pOkpySdWurbWtK/pP+2n5D0I0m1pbitJV0raaOkJ/La+t22SnwlXf8lkk4a6nIdJgchqQL4GvBaYB7wdknzsq2qaLqASyNiHnAK8OF0XS8DfhMRc4HfpO9LzT8BT+W9/xzwpYg4AtgCvD+TqorrP4A7IuJo4HiS9S/ZbS1pOvCPwPyIOBaoAC6kNLf1d4Fz+rQdaNu+FpibPi4GrhrqQh0mB7cAWB4RKyKiA7gBOC/jmooiItZFxCPp6+0kv1ymk6zvdels1wFvzKTAIpE0A/gb4NvpewFnAjels5TiOjcCrwCuAYiIjohoo8S3NcltyuskVQL1wDpKcFtHxO+BF/o0H2jbngdcH4kHgCZJU4eyXIfJwU0HVue9X5O2lTRJs4ATgQeBKRGxLp20HpiSVV1F8mXgY0BP+v4QoC0iutL3pbjNZwOtwHfS7r1vSxpHCW/riFgL/DvwHEmIbAUWUfrbuteBtu2w/Y5zmNh+JI0Hbgb+OSK25U+L5DjykjmWXNK5wMaIWJR1LSOsEjgJuCoiTgR20qdLqwS39USSv8JnA9OAcby4K6gsFGvbOkwObi0wM+/9jLStJEmqIgmSH0TET9PmDb27venzxqzqK4K/At4gaSVJF+aZJGMJTWlXCJTmNl8DrImIB9P3N5GESylv67OAZyOiNSI6gZ+SbP9S39a9DrRth+13nMPk4B4G5qZHfFSTDNjdmnFNRZGOFVwDPBUR/y9v0q3ARenri4Cfj3RtxRIRl0fEjIiYRbJtfxsRfwvcDZyfzlZS6wwQEeuB1ZKOSpteBTxJCW9rku6tUyTVp//We9e5pLd1ngNt21uBd6dHdZ0CbM3rDhsUnwH/F0h6HUm/egVwbURckW1FxSHpdOBe4HH2jR98gmTc5EbgMJJL9r8tIvoO7o15kl4JfDQizpU0h2RPZRLwKPDOiGjPsLxhJ+kEkoMOqoEVwHtJ/rgs2W0t6bPABSRHLj4K/B3J+EBJbWtJPwJeSXKp+Q3Ap4Gf0c+2TYP1qyRdfruA90bEwiEt12FiZmaFcjeXmZkVzGFiZmYFc5iYmVnBHCZmZlYwh4mZmRXMYWJWIEl/TJ9nSXrHMH/3J/pbltlo40ODzYZJ/rkqg/hMZd61ofqbviMixg9DeWZF5T0TswJJ2pG+vBI4Q9Jj6b0zKiR9QdLD6b0i/j6d/5WS7pV0K8lZ2Ej6maRF6f02Lk7briS5yu1jkn6Qv6z0jOUvpPfmeFzSBXnf/bu8e5X8ID0xzayoKv/yLGY2QJeRt2eShsLWiDhZUg3wB0m/Suc9CTg2Ip5N378vPSO5DnhY0s0RcZmkSyLihH6W9WbgBJJ7kTSnn/l9Ou1E4BjgeeAPJNegum+4V9Ysn/dMzIrnNSTXPXqM5LI0h5DchAjgobwgAfhHSYuBB0guvDeXgzsd+FFEdEfEBuAe4OS8714TET3AY8CsYVgXs4PynolZ8Qj4h4i4c7/GZGxlZ5/3ZwGnRsQuSb8DagtYbv61pbrx/3MbAd4zMRs+24GGvPd3Ah9KL+2PpCPTm1D11QhsSYPkaJLbJvfq7P18H/cCF6TjMi0kd058aFjWwmwI/BeL2fBZAnSn3VXfJbk3yizgkXQQvJX+bwt7B/BBSU8By0i6unpdDSyR9Eh6efxetwCnAotJbnT0sYhYn4aR2YjzocFmZlYwd3OZmVnBHCZmZlYwh4mZmRXMYWJmZgVzmJiZWcEcJmZmVjCHiZmZFez/A8/ym04GWewWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 记录优化中间过程\n",
    "loss_list, singular_value_list = [], []\n",
    "U_learned, V_dagger_learned = [], []\n",
    "    \n",
    "# 确定网络的参数维度\n",
    "net = NET(M, weight)\n",
    "\n",
    "# 一般来说，我们利用 Adam 优化器来获得相对好的收敛\n",
    "# 当然你可以改成 SGD 或者是 RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# 优化循环\n",
    "for itr in range(ITR):\n",
    "\n",
    "    #  前向传播计算损失函数\n",
    "    U, V_dagger, loss, singular_values = net()\n",
    "\n",
    "    # 反向传播极小化损失函数\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # 记录优化中间结果\n",
    "    loss_list.append(loss[0][0].numpy())\n",
    "    singular_value_list.append(singular_values)\n",
    "    \n",
    "    if itr% 10 == 0:\n",
    "        print('iter:', itr,'loss:','%.4f'% loss.numpy()[0])\n",
    "        \n",
    "# 绘制学习曲线\n",
    "loss_plot(loss_list)\n",
    "\n",
    "# 记录最后学出的两个酉矩阵    \n",
    "U_learned = U.numpy()\n",
    "V_dagger_learned = V_dagger.numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接着我们来探究下量子版本的奇异值分解的精度问题。在上述部分，我们提到过可以用分解得到的更少的信息来表达原矩阵。具体来说，就是用前 $T$ 个奇异值和前 $T$ 列左右奇异向量重构一个矩阵：\n",
    "\n",
    "$$\n",
    "M_{re}^{(T)} = UDV^{\\dagger}, \\tag{4}\n",
    "$$\n",
    "\n",
    "并且对于一个本身秩 (rank) 为 $r$ 的矩阵 $M$, 误差随着使用奇异值的数量变多会越来越小。经典的奇异值算法可以保证：\n",
    "\n",
    "$$\n",
    "\\lim_{T\\rightarrow r} ||M - M_{re}^{(T)}||^2_2 = 0, \\tag{5}\n",
    "$$\n",
    "\n",
    "其中矩阵间的距离测量由 Frobenius-norm 来计算,\n",
    "\n",
    "$$\n",
    "||M||_2 = \\sqrt{\\sum_{i,j} |M_{ij}|^2}. \\tag{6}\n",
    "$$\n",
    "\n",
    "目前量子版本的奇异值分解还需要很长时间的优化，理论上只能保证上述误差不断减小。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:49:24.407257Z",
     "start_time": "2021-03-09T03:49:24.027591Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABFGElEQVR4nO3dd3hUZfbA8e9JIYFQAiGEJr3XBKKADUGaK6KLFFlFUVfXviuLbf2t4rquHV1dV9eyFlSKCCjqinRUBKV3pAWkhhYkEEg7vz/emxBCygQmmYScz/PMk5l779x7ZpLMmbeLqmKMMab8Cgp0AMYYYwLLEoExxpRzlgiMMaacs0RgjDHlnCUCY4wp50ICHcCZqFmzpjZq1CjQYRhjTJmyZMmS/aoanXt7mUwEjRo1YvHixYEOwxhjyhQR2ZbXdqsaMsaYcs4SgTHGlHOWCIwxppwrk20ExpQGiYmJjBo1ivXr15OZmRnocIwBICgoiFatWvHCCy9Qq1Ytn55jicCYMzRq1Ch69OjBO++8Q2hoaKDDMQaAtLQ0xo4dy6hRo/jggw98ek65qRqaumwnA56exKLHunDV058yddnOQIdkyrj169dzww03WBIwpUpoaCjDhw9n/fr1Pj+nXJQIpi7bySOTV/EXHcf5wRsYfPRjHpkcAcA1cfUCHJ0pqzIzMy0JmFIpNDS0SNWV5aJE8Pz0DVRO28/g4HkEiTI4eD6V0w7w/PQNgQ7NGGMCrlwkgl1JKdwXMoVgMgAIIoN7QyazKyklwJEZc3aCg4OJjY2lXbt2XHXVVSQlJQUslrlz57JgwQK/nW/q1KmsXbs2+/Fjjz3GzJkz/Xb+3D7//HOeeeYZn449duwYUVFR/Prrr6dsv+aaa5gwYQLg4u/QoQOtWrWiXbt2TJo0Kfu4hQsX0qVLF2JjY2ndujWjR48mISGB+vXrn/ZNPjY2lkWLFjF69Gjq1atHbGwszZs3Z+DAgae8P2ejXCSC9tVSGBw8j1Bxb3CYZDA4eD7tqh0PcGSmPJm6bCcXPTObxg9/yUXPzPZLO1XFihVZvnw5q1evpkaNGrz22mt+iPTMFJQI0tPTi3y+3Ingb3/7G7169Trj+AozYMAAHn74YZ+OrVSpEn379mXKlCnZ2w4fPsx3333HVVddxYoVKxg1ahSfffYZ69evZ9q0aTz00EMsWbIEgJtuuok333wz+3c3ZMgQGjVqRIMGDfj222+zz7l+/XqOHDlCly5dALj//vtZvnw5GzduZOjQofTs2ZN9+/ad9WsvF4ng5TozEE5diS2ITP5Z55sARWTKm6x2qp1JKSiwMymFRyav8munhW7durFzpzvf5s2b6devH507d+aSSy7Jbjjcu3cvv/3tb+nYsSMdO3bM/uAeM2YM7dq1o127drz88ssAJCQk0Lp1a2677Tbatm1Lnz59SElxpehXXnmFNm3a0KFDB6677joSEhJ44403eOmll4iNjeXbb79lxIgR3HHHHXTp0oUHH3yQ0aNH88ILL2TH265dOxISEgD44IMP6NChAx07dmT48OEsWLCAzz//nAceeIDY2Fg2b97MiBEjsr9Vz5o1i7i4ONq3b88tt9zCiRMnADf9zOOPP06nTp1o3759ng2mXbt2Zc2aNdmPL7vsMhYvXsx7773HPffcA8C0adPo0qULcXFx9OrVi7179552nmHDhjF+/Pjsx1OmTKFv375UqlSJF154gb/85S80btwYgMaNG/OXv/yFF198EXBdj+vUqQO4Ul2bNm3yPOf48eO57rrr8vx9Dx06lD59+vDxxx/nub8oykVjcZPja0BO/UYSJuluuzF+MvQ/P+S7b9n2JFIzTi3yp6RlMHraGq6Jq8fBo6nc+eGSU/ZP+EM3n6+dkZHBrFmzuPXWWwG4/fbbeeONN2jevDmLFi3irrvuYvbs2dx33310796dKVOmkJGRQXJyMkuWLOHdd99l0aJFqCpdunShe/fuVK9enY0bNzJu3DjeeusthgwZwqeffsoNN9zAM888w9atWwkLCyMpKYnIyEjuuOMOKleuzKhRowB455132LFjBwsWLCA4OJjRo0fnGfuaNWv4+9//zoIFC6hZsyYHDx6kRo0aDBgwgP79+zNo0KBTjj9+/DgjRoxg1qxZtGjRghtvvJHXX3+dP/3pTwDUrFmTpUuX8u9//5sXXniBt99++5TnDx06lIkTJ/LEE0+we/dudu/eTXx8PKtXr84+5uKLL2bhwoWICG+//TbPPfdc9od4lr59+/L73/+eAwcOEBUVxfjx47MTyZo1a7Lfhyzx8fG8+uqrgPtm37JlSy677DL69evHTTfdRHh4OEOGDCE2NpZXX32VkJAQJkyYwCeffJLv771Tp05F6h2Un3JRIuCO72D0YRh9mKR7NnBCQ1lYc6DbbkwJyJ0EsiQdSzur86akpBAbG0vt2rXZu3cvvXv3Jjk5mQULFjB48GBiY2P5wx/+wO7duwGYPXs2d955J+C+iVarVo3vvvuO3/72t0RERFC5cmUGDhyYXT3RuHFjYmNjAejcuXP2N/gOHTpw/fXX8+GHHxISkv/3ycGDBxMcHFzga5g9ezaDBw+mZs2aANSoUaPA4zds2EDjxo1p0aIF4KpZ5s+fn71/4MCBp8Wb05AhQ7JLFhMnTjwt0QDs2LGDvn370r59e55//vlTShBZKlSowIABA5g0aRL79+9n2bJl9O3bt8DYszz22GMsXrw4+xt9v379AIiJiaFdu3bMmjWL5cuXExISQrt27fI9j7/WnC8XJYKcImvW5pM69zF1d3Xan0gnIqzcvQWmmBT0Df6iZ2azM4/OCfUiKwJQI6JCkUoAWbLaCI4dO0bfvn157bXXGDFiBJGRkSxfvrzI58stLCws+35wcHB21dCXX37J/PnzmTZtGk899RSrVq3K8/kRERHZ90NCQk5pCD1+vHja6LJiDg4OzrNtol69ekRFRbFy5UomTJjAG2+8cdox9957LyNHjmTAgAHMnTs339LMsGHDePLJJ1FVrr766uzuxG3atGHJkiV07Ngx+9glS5YQHx+f/bhp06bceeed3HbbbURHR2eXLLKqh2JiYhg2bFiBr3XZsmWnnPNMlWiJQEQiRWSSiKwXkXUi0k1EaojIDBHZ6P2sXtxxNLniXr4/0YTPV+wq7ksZA8ADfVtSMfTUb8YVQ4N5oG9Lv5y/UqVKvPLKK7z44otUqlSJxo0bZ1cpqCorVqwA4PLLL+f1118HXHXS4cOHueSSS5g6dSrHjh3j6NGjTJkyhUsuuSTfa2VmZvLLL7/Qo0cPnn32WQ4fPkxycjJVqlThyJEj+T6vUaNGLF26FIClS5eydetWAHr27Mknn3zCgQMHADh48CBAvudr2bIlCQkJbNq0CYCxY8fSvXv3Ir1fQ4cO5bnnnuPw4cN06NDhtP2HDx+mXj03xuj999/P9zyXXXYZGzdu5LXXXjvlQ3vUqFE8/fTT2SWShIQEXn75ZR544AHAJdKsb/MbN24kODiYyMhIwJVovvrqKyZMmJBv+wDAp59+yjfffFNosvBFSVcN/RP4WlVbAR2BdcDDwCxVbQ7M8h4Xq04NIukdncSJOc/7rWhlTEGuiavH0wPbUy+yIoIrCTw9sL1fBzTGxcXRoUMHxo0bx0cffcQ777xDx44dadu2LZ999hkA//znP5kzZw7t27enc+fOrF27lk6dOjFixAguuOACunTpwu9//3vi4uLyvU5GRgY33HAD7du3Jy4ujvvuu4/IyEiuuuoqpkyZkt1YnNu1117LwYMHadu2Lf/617+yq3batm3Lo48+Svfu3enYsSMjR44E4LrrruP5558nLi6OzZs3Z58nPDycd999l8GDB9O+fXuCgoK44447ivReDRo0iPHjxzNkyJA8948ePZrBgwfTuXPn7CqrvAQFBTFo0CAOHDhwSjKKjY3l2Wef5aqrrqJFixa0aNGC119/nZYtXeIfO3YsLVu2JDY2luHDh/PRRx9lV6FFRkbSrVs3YmJiaNKkySnXy2qMb968OR9++CGzZ88mOvq0dWaKTErqg1BEqgHLgSaa46IisgG4TFV3i0gdYK6qFvg1KT4+Xs92YZpFE5+ly9p/sHXgFzTukP+3H2PyEx8fbwskGZ88/PDDLFq0iOnTp1OhQoUSuWZef58iskRVT6tLKskSQWNgH/CuiCwTkbdFJAKIUdXd3jF7gJi8niwit4vIYhFZ7I9+s2373kZmSEUab51w1ucyxpiCPPPMM8yZM6fEkkBRlWQiCAE6Aa+rahxwlFzVQF5JIc8iiqq+qarxqhrvj6JQ5Wo1COowGFZ/CilJZ30+Y4wpq0oyEewAdqjqIu/xJFxi2OtVCeH9TCypgFLjRkDaMRZPO73XgDHGlBcllghUdQ/wi4hk1f9fDqwFPgdu8rbdBHxWUjFVOK8zayp0IP3Yr4UfbIwx56iS7kR/L/CRiFQAtgA345LRRBG5FdgG5N2MX0zaPDwPCSof4+qMMSYvJZoIVHU5kNfoh8tLMo6cJCgIzcxk59Z11G/aNlBhGGNMwNhXYeDHsY8S/UF3DiTaADNTttg01P5TlGmowU1Fff3119O+fXvatWvHxRdfTHJyMj169GD69OmnHPvyyy9z5513kpCQQMWKFYmLi6N169ZccMEFvPfee35+JUVniQCoc8FAwiSNDdPfDHQo5lx3ZA+8ewUcOX02yzNh01D7T1GmoQY3OC8mJoZVq1axevXq7LWrc88gCm4W0awRwE2bNmXZsmWsW7eO8ePH8/LLL/Puu+/69bUUlSUCoEHr81kX2pbztkwgMyMj0OGYc9m852D7Qpj3rN9PbdNQl+w01Lt3786ehgLc1BdhYWEMGjSIL7/8ktTU1Oz3cdeuXXlO29GkSRPGjBnDK6+8Utivt3ipapm7de7cWf1t8bQ3VB+vqivmTfX7uc256bS/w//+5vTbojfdvhNHVd/spTo6UvXxqu7nW71Ul37o9ifvP/25PoiIiFBV1fT0dB00aJD+73//U1XVnj176s8//6yqqgsXLtQePXqoquqQIUP0pZdeyn5OUlKSLl68WNu1a6fJycl65MgRbdOmjS5dulS3bt2qwcHBumzZMlVVHTx4sI4dO1ZVVevUqaPHjx9XVdVDhw6pqurjjz+uzz//fHZsN910k1555ZWanp6e5/62bdvq1q1bdfXq1dq8eXPdt2+fqqoeOHAg+/mffPLJKef75JNPNCUlRevXr68bNmxQVdXhw4dnv6aGDRvqK6+8oqqqr732mt56662nvWdjxozRxx57TFVVd+3apS1atFBV1XfffVfvvvtuVVU9ePCgZmZmqqrqW2+9pSNHjjztPMuWLdPo6Gjt2rWrPvroo9nvt6rqlVdeqVOnus+Sp59+Wv/85z+rqurWrVu1bdu2p5zn0KFDGh4eftr5z1Zen5PAYs3jM9VKBJ52vW4giSoc/XFsoEMx56rD2yFrdhVVSNp+1qe0aagDNw11bGwsW7Zs4YEHHuDgwYOcf/75rFu3Djh1gZmc1UJ50VIw35nNwewJC4/g43av8vxSYfbh49SuFh7okExZc/OX+e878SscT+LkwHl1j5t5dd4RUQU/Px82DXX+MZfENNRZiXPgwIEEBQXx1Vdf0bp1a66++mruv/9+li5dyrFjx+jcuXO+8S5btozWrVuf2Yv1EysR5HB5j74cywxl/E9n/03NmFPMew401+I0mum3tgKbhrrkp6H+/vvvOXToEACpqamsXbuWhg0bAi5B9OjRg1tuuaXA0kBCQgKjRo3i3nvvLVL8/maJIIcGUZUYWX89F31/C+lpZ7dylDGn2PEjZKSeui0j1W33E5uG2nf+mIZ68+bNdO/ePft9iI+P59prr83eP2zYMFasWHFaIti8eXN299EhQ4Zw3333cfPNNxcpfn8rsWmo/ckf01DnZ/n094j94Y/s+s371L3gmmK5hjk32DTUpjQrrdNQlwnte/4OrRxD3U3jAh2KMcaUCEsEuQSHVkDihqM/T+dY4tZAh2OMMcXOEkEejnccjgLLpv4z0KGYUiwoKIg0a0sypVBaWhpBRZhM0xJBHsJrNuKnBr+ncvOLAx2KKcVatWrF2LFjLRmYUiUtLY2xY8fSqlUrn59jjcXGnKHExERGjRrF+vXrT+kfb0wgBQUF0apVK1544QVq1ap1yr78GottQFkBEnduZeuPX9Llt/cEOhRTCtWqVYsPPvgg0GEYc9asaqgACTP+Q5cVj7Jt48pAh2KMMcXGEkEBmvW7k3QNYsdMW9PYGHPuskRQgBq1G7K6ysW02fs5KceOBTocY4wpFpYIChHW9fdU5wgrZlhdsDHm3GSJoBCtLuzPtqD6bNuwPNChGGNMsbBEUAgJCmZuj8k8dPAqVu88HOhwjDHG7ywR+OCazk0IDw3i0wWnL05hjDFlnSUCH1SrFMqYOrP50+pBHPk1KdDhGGOMX1ki8FGL83tTTY6SserTQIdijDF+ZYnAR83ie0N0KyLXfBjoUIwxxq9KNBGISIKIrBKR5SKy2NtWQ0RmiMhG72f1kozJZyIQfwvsWsrOtT8EOhpjjPGbQJQIeqhqbI6Jjx4GZqlqc2CW97hUOt5mMClage3fvBboUIwxxm9Kw6RzVwOXefffB+YCDwUqmIKEV6nB6p5v0KrNhYEOxRhj/KakSwQKfCMiS0Tkdm9bjKru9u7vAWLyeqKI3C4ii0Vk8b59+0oi1jy1634t1aPrBOz6xhjjbyWdCC5W1U7AFcDdInJpzp3qFkfIc4EEVX1TVeNVNT46OroEQs3fqrmTWDhmKGpz0BtjzgElmghUdaf3MxGYAlwA7BWROgDez8SSjOlMpB3cTtdfv2bVjzMDHYoxxpy1EksEIhIhIlWy7gN9gNXA58BN3mE3AZ+VVExnqk3fW0mmIse+fyvQoRhjzFkrycbiGGCKiGRd92NV/VpEfgImisitwDZgSAnGdEbCI6qxOPoK4hKnsS9xN9G1rM3AGFN2lViJQFW3qGpH79ZWVZ/yth9Q1ctVtbmq9lLVgyUV09mI6XknYZLGhulvBjoUY4w5K0VKBCISLyJDvaqdrOqe0tAFtcSd1/oC5lTqx9fbg8jIzLN92xhjygSfEoGIxIjIQuBH4GNOdvEcA7xYTLGVesf6vcyHRzox/+fAdWc1xpiz5WuJ4CVgLxAF5Fyz8RNco2+51LtNDA0jMlg+Z1KgQzHGmDPma7XO5cDlqnrIa+zNshlo4PeoyogKIUE8X+tr4naNZ/fO/tSpV27fCmNMGeZriaAikJrH9mjguP/CKXsa9r6TUMmg5saJgQ7FGGPOiK+JYD4wIsdjFZFg3JxAs/wdVFkS06Q9NLqE0OXvg400NsaUQb4mggeB20RkBhCGayBeC1wEPFJMsZUZaZ1GQNJ2Vn87OdChGGNMkfmUCFR1LdAeWAB8A4TjGorjVHVz8YVXNgS1voqDVOPw2rmBDsUYY4rM5zEAqroHeLwYYymzgkPDCL5nERfVtBHGxpiyx9dxBPeIyA15bL9BRO7yf1hlTzUvCaSeOBHgSIwxpmh8bSP4E/BLHtsTgPv9FUxZ99PHT5D4dHuOWzIwxpQhviaC+rgJ4XLb4e0zQOU6LajPXlbMGh/oUIwxxme+JoI9QGwe2zsB+/0WTRnX8pJB7JMowpa/H+hQjDHGZ74mgo+BV0Skt4iEerc+wMvAR8UWXRkTFBLK9kaDiU1dwpafVwc6HGOM8YmvieBx4HtgOm6uoWPA/3DdSf9aPKGVTU363kW6BrFz5uuBDsUYY3ziU/dRVU0DhonIY0Acbl3h5aq6sTiDK4uq127Ix3Uf5MNddemcmk6lCuVylm5jTBlSpPUIVHWjqk5U1U8sCeSvRd8/sPZETaat2BXoUIwxplA+f10VkaG4WUhrkSuBqOoAP8dVpnVuWJ1ra26H2Z/D+f8JdDjGGFMgXweUPQ98CDQCkoADuW4mBxHh+jq7GZoynp0/Lw90OMYYUyBfSwQ3AsNU1VZg8VHLK+5EN/+belsmQIvYQIdjjDH58rWNIAhYXoxxnHMiatRB2gyA5R+hqccKf4IxxgSIr4ngTeC0uYZMwY53vBGOH+b7z98OdCjGGJMvX6uGIoHfiUhvYCWQlnOnqt7n57jOCeHNurO2YidCgqTwg40xJkB8TQRtOFk11CrXPvVbNOcaEdo8NCfQURhjTIF8HVDWw18X9Ja4XAzsVNX+ItIYGA9EAUuA4aqa1/rIZVZa6gm2rF1My9iLAh2KMcacpkgDyvzkj8C6HI+fBV5S1WbAIeDWAMRUrFa+czf1pgzkwEHraWuMKX18TgQi0kNE3hSRr0Vkds5bEc5RH7gSeNt7LEBPIKtb6vvANT5HX0ZEX3gDleU4a762RmNjTOnj64CyEbhJ5qoAlwH7gOq4aajXFuF6LwMPApne4yggSVXTvcc7gHr5xHC7iCwWkcX79u0rwiUDr0GH7iSENKHOpo/JzMgs/AnGGFOCfC0RjALuUdVhuB5Dj6hqHG60cbIvJxCR/kCiqi45k0BV9U1VjVfV+Ojo6DM5ReCIcLjtDTTPTGDlj7MCHY0xxpzC10TQBJjp3T8BVPbu/wsY4eM5LgIGiEgCrnG4J/BPIFJEshqt6wM7fTxfmdKqz60cJZzEhRMDHYoxxpzC10RwAFctBO6Dup13Pwqo6MsJVPURVa2vqo2A64DZqno9MAcY5B12E/CZjzGVKWERkXzY4X3u2nc1ew4fD3Q4xhiTzddE8C3Qx7s/Ebda2bvAOGDGWcbwEDBSRDbhEss7Z3m+Uqtf90tIzxQm/Lg90KEYY0w2XweU3QOEe/efBtJxVT0Tgb8X9aKqOheY693fAlxQ1HOURQ2jInii9gJiF/yN9B4LCAkJDnRIxhjj84CygznuZ+L6/psz0KlZPdovXk/iujnUat8r0OEYY4zP3UczRKRWHtujRCTD/2Gdu9r2vgkNr0atDeMCHYoxxgC+txHkN2taGHBOTQdR3ILCIpCOw9C1n/HrPlvK0hgTeAVWDYnISO+uAneISM4xA8HAJcD6YortnHW0w3AiFr3B0mmvcdktTwU6HGNMOVdYG8G93k8Bfg/krAZKBRKAO/wf1rktol47Fje9h3rt+wU6FGOMKTgRqGpjABGZAwxU1UMlElU5ED/cSgLGmNLBpzYCVe2ROwmISDMRCc/vOaZwW1cv5LuPnwl0GMaYcs7XXkP/EJGbvPsiIjOBn4HdItKlOAM8lx1a9DFdNzzL1i0bAx2KMaYc87XX0PXABu/+FUBHoCvwAWBfac9Qwz53ESKZbJv5RqBDMcaUY74mghjcFNEAvwEmquqPwKtAXHEEVh5EndeKdZXiabVrCsdPnAh0OMaYcqook8419O73AbLmUg4h/zEGxgdB599KbQ6wdOaEQIdijCmnfE0EnwIfi8gMoAYw3dseC2wqhrjKjRaXDGJbUH1WrrfhGMaYwPA1EYwEXsGtRtZbVY962+sArxdHYOWFhFRgZo/PeWbfRazZdTjQ4RhjyiFfu4+mq+qLqvpHVV2WY/tLqmoL8Z6lazufR1iI8OW3PwU6FGNMOZTvgDIR6QQsV9VM736+VHWp3yMrRyIrVeCd6E9otXYWycfWU7lSpUCHZIwpRwoaWbwYqA0keveVvBuGFTfvkDkL9S4YQM3pk0le/yV0GhzocIwx5UhBiaAxsC/HfVOMGncZAAsbUHnVB5YIjDElKt9EoKrb8rpviklQMJmdbiJozpMkrF9Oo1axgY7IGFNO+DrFRBMRGSki/xKRV0XkfhGxUoKfpbQbRpoGs232W4EOxRhTjhS6VKWI/Bm3TnEwrr1AgGjgWRF5SFVfKt4Qy4+IqHps6j+eru0vDnQoxphypMASgYhcDDwHPA9Eq2odVa0N1AJeBJ4XkYuKP8zyo9n5fQgLt15DxpiSU1jV0J3AB6r6aK4F7A+o6iPAh8BdxRlgebTss1f54blrUNVAh2KMKQcKSwRdgfcK2P+ed4zxo+ATh+l2bA53PPkSix7rwlVPf8rUZTsDHZYx5hxVWCKoDWwpYP9m3DQTxo+21r+aExrK/envcL5sYPDRj3lk8ipLBsaYYlFYIqgIFDQ/cioQ5suFRCRcRH4UkRUiskZEnvC2NxaRRSKySUQmiEgF30I/dz03fx8zM+NoKTsIEmVw8Hwqpx3g+ekbCn+yMcYUUaG9hoArRSS/2dAii3CtE0BPVU0WkVDgOxH5H25Cu5dUdbyIvAHcSjmfyG5XUgoScrJ9IIhM7g2ZzONJtwQwKmPMucqXRPBOIft9atFU1/KZ7D0M9W4K9AR+521/HxhNOU8E7aul0PP4csSb0CNM0hkSPI8DQTXYvKsLTevWDGyAxphzSoFVQ6oa5MPN53mGRCRYRJbjxiPMwLUxJKlqunfIDqBePs+9XUQWi8jiffv25XXIOePlOjOQXPk1hAzuD5pIzHvdYPF/ISMtQNEZY841vq5H4BeqmqGqsUB94AKgVRGe+6aqxqtqfHR0dHGFWCo0Ob6GMEk/ZVuIZJIe2YiK0Q3hi/s5NiaO8W89S/JxSwjGmLPjS9WQ36lqkojMAboBkSIS4pUK6gPWNeaO7/LcHAKgChtncOzzv9J67xdUrPAgACfSMwgLsUlgjTFFV2IlAhGJFpFI735FoDewDpgDDPIOuwn4rKRiKpNEoEUfao78gXZ/mkpwkJCcmMDPT3Xl47Fvkng4JdARGmPKmJKsGqoDzBGRlcBPwAxV/QJ4CBgpIpuAKApvnDYAQUEEV44CICNpJ3VDkvnd5gfYNeYSPhr3PgeSC+r1a4wxJ0lZnMYgPj5eFy9eHOgwSpeMNPZ/9y7B3z5H9fR9LND2fN/1DX5/aQuqR5T7oRnGGEBElqhqfO7tJdpYbIpRcCg1u99O9YdWs+/iv3GsRlv+PX8blzw3h3emzeZwijUqG2Py5lNjsYhUx/Xv74GbefSUBKKqtfwemTkzoeFE9/ojvXrB13uOMPmLL7h58Z0sXXcx8SNegFo+d9QyxpQTvvYa+gBoixvwtRcfB5GZwGpZuwqP3PAbEr+5j7hVb8O/u3K89bVMrjaca3peTKUKAek0ZowpZXxqIxCRI0B3VV1a/CEVztoIzsDRA/D9y6QvfJMjGSEcuWsFDWJshLIx5cnZthFsLsKxpjSKiII+TxJy/wrSBrzhkoAq/3vjISbOXcKJ9IxAR2iMCRBfP9z/CDwtIh1FxEYtlWVValOr81UAnNixnD573qT/nCuY8PTvmfTtKtIyMgMcoDGmpPmaCDbhpqReCqSKSEbOW/GFZ4pT2HlxBN3zI8mN+nBDxhT6zOzNB8/cyeRFP5NuCcGYcsPXNoL5QHXgDfJoLFbVT4slunxYG4H/6Z7V7J/2OJm7V3LJseepXzOSP17ejP4d6xEcJIEOzxjjB/m1EfjabSQeuEBVV/s3LFNaSO12RN/2KZqSxKtbjvOvb1bTcMpVvD2rL7fdN5qgCuGBDtEYU0x8rRpaC1QtzkBM6SAVI+nbtjafjWjJebWi+EPy6wS9Fo8ueZ/vf95DWRyJbowpmK+J4P+AMSLSS0RiRKRGzltxBmgCI6j6eUTdPQOGT4HKtZBp91Hnw+7M/HFFoEMzxviZr1VDX3k/v+HU9gHxHltPonORCDTtCU16kLHuS4IXTuKyzu0B+HbB91CzORc3j+az5bt4fvoGdiWlUDeyIg/0bck1cXmuL2SMKYV8TQQ9ijUKU7qJENymPw3b9AdAk/cRP2MgGzPq8HClG5mS3Ipq6QcZX+FV7km6j0cmpwJYMjCmjCi011DWQvPAjaq6oUSiKoT1GgqwjHTSlo/j+Mx/UCVlF4syW5GUGUHv4KV8mHE5j6XfQr3Iinz/cM9AR2qMyeGMRxarahrQGJtfyGQJDiG083Cq/HkFf027mSayi74hSwgSZXDwfKJJYmdSCpsSjwQ6UmOMD3xtLH4fuK04AzFlUEgFZlcZwIyMzqSpayYKIpMPKzzFAyHjSd2/FYB5P+/jj+OX2WI5xpRSvrYRRADXi0hvYAlwNOdOVb3P34GZsuH/utegx9ffEeoNMA+TdJqxixYhu2DiNGjWi5DqA1i1vS5VwkMBeHnmz6zccZj4RtU5v1EN2terRnio9TcwJlB8TQStcdNLADTJtc+qjMqxKw58QEYQp/4VBIUg7a+F6o1g6ftctGkGsy9/HEJ6ARAWEsz2g8eYvT4RgArBQXQ8rxrxjWpwfqPqdG5Qg2qVQkv8tRhTXtlSlebsvHEx7Fl1+vba7eGO7yAjHX7+Gup1gqp1YcP/YNmHcP6tHKjVjSXbD7N42yF+SjjIqh2HSc90f4+928Tw1o2uTevQ0VRbbtMYPzjbKSayThIONMN9/9usqsf9FJ8pq+74ruD9wSHQuv/Jx8cOwvYfYP0XRFVvTJ/4m+nT/Xr4TWtSUjNYsSOJxQkHsxfNycxUerw4l2s71eev/dugqmzYe4QWtaoQZHMgGeMXvk46Fwr8A7gHqIAbSHYCeBV41OtZVGKsRFDGpZ+AddNg8X9h2/dQsyXcvcgNYMslNT2TcT9up0VMFbo1jWJTYjK9xsyjangInRtW96qTatChvrUzGFOY/EoEviaCMcAw4GHcmAKAS4CngY9UdZQfYy2UJYJzSOI6OLIHmvaAtOMw9rfQ9hroMBQqRp52+OGUNGat28tPCYdYnHCQjYnJgGtn6FD/ZDtDlyZRVA6zpTiNyelsE8Ee4BZV/SrX9iuBt1W1jt8i9YElgnPUwS0w6VbYtRRCK0G7ayH+Fte+kN9TjqayZJtLCj8lHGTVzsOkZSif3NGN8xvVYO2uX/l57xH6tatNeGgwU5fttOkwTLl1tm0E1XDLVea2GYg8i7iMOalGE7h9DuxaBovfhVWfwLKxcPtcqBuX91MiKtC7TQy928QAcDwtg+W/JNG+XjUAvli5i7e+3UK/drWZumwnD0xaQVqG+/KzMymFRya7hm5LBqY887VEsBBYoqp359r+OhCrqt18OMd5wAdADK6x+U1V/ac3e+kEoBGQAAxR1UMFnctKBOXE8cOul1GHoa79YMZjkHoM4m+GmLY+nSIjU9l+8BiNa0Zw0TOz2ZmUctoxdaqF88Mjl/s7emNKnbNdvP5B4CYR2SAi73u3DcANwAM+niMd+LOqtgG6AneLSBtcu8MsVW0OzPIeGwPh1aDjdScbkU8kw9IP4PUL4Z2+sGKCa1coQHCQ0LhmBAC78kgCALsPH+dP45excMsBW2/BlEs+JQJVnQ+0ACYBlb3bJ0BLVS2k/2D2OXar6lLv/hFgHVAPuBo3hQXez2uKEL8pT/qPgZHroM/f4WgiTLkdZo72+el1IyvmuT2iQjCz1iVy3ZsLufzFecxYu9dPARtTNvjcrUJVdwGP+uOiItIIiAMWATGqutvbtQdXdZTXc24Hbgdo0KCBP8IwZVFEFFx4L3S9GxLmQ7Xz3Pbti2DuP1zjcsvfQPDpI5Mf6NuSRyavIiUtI3tbxdBgnvpte/q0jeHLlbsZ/9MvhAa7EsjOpBS27EvmoqY1bcyCOacV2Ebg6+pjqnrQ5wuKVAbmAU+p6mQRSVLVyBz7D6lq9YLOYW0E5jTrv4KvHoBfd0DlGOh0I3S6CSLPO+WwovQaevGbDfx77mZ+eLgntaqGk5GpBFtCMGXYGXUfFZFMCp9LSFXVp5KFNzDtC2C6qo7xtm0ALlPV3SJSB5irqi0LOo8lApOnzAzYOMMNVNv4DUTUhJHr3ejmM3AiPYPl25Po0iQKgOHvLCIsJJhhF5xH9xbRhAT72sRmTOlwpt1HC1qZrB/wR1wjsC8BCPAOsC4rCXg+B24CnvF+fubL+Yw5TVAwtOznbknbYd8GlwQyM2HsNdD4UogbDihMuhkGvQdV8qyJBNzkeFlJIDNTaV+vGhMX72Dmur3UrhrOkPj6DDn/POpXr1QiL8+Y4lLkSedEJA54Hjey+D/Ak6q6z4fnXQx8C6wCMr3Nf8G1E0wEGgDbcN1HC6xqshKBKZKj+2HSLbB1HgSFQLX6cGiba0/oP6bw5+eQlpHJrHWJjP9pO/N+dn/2lzSPZtj553F56xgqhFgpwZReZzWy2DtBY+ApYDAwGfiLquY1yKzYWSIwZ2T/JvjhX7DkXfc4uAL8aTWEhEFIOISGF+l0O5NSmPjTL0xc/Au7Dx+nZuUKvHVjPHENCmziMiZgznhksYhEAY8BdwDfAxeq6k/+D9GYYlazGUiQSwAZqW7bvGehYnVY9AY06wWt+kOLPm4MQyHqRVbk/t4tuO/y5sz/eR+Tlu6gWa3KAMxcu5ejqelc1aGu9TgypV6BiUBEHsUNGEsArlbVr0siKGOKxZE9sPyjk0kgI9U9HvQutB8MG76CtVMhKBRa9IWhH+Y5I2puwUFCj1a16NGqVva28T/9wu7DKVwd63okJR45Tq0qRStxGFNSCisRPAmkADuAu0TkrrwOUtUB/g7MGL+b9xxo5qnbNBM2zYSrXoYrx8COn2D9F5CRdjIJTL0LarZwpYWazXy61JvDO7PfW6P5yPE0uj83l9Z1qnDdBQ3o36FO9noLxpQGhf01foAtRWnOFTt+PFkayJKR6rYDBAVBgy7uliX1GCSudSWHmY9DdCtodSV0uA6iW+R7qaAgoVZVVwIQEf7cpwXjftzOg5NW8uS0tQyIrcuwCxrQrl7hVVDGFDdbqtIYXyT94qqO1n8BCd/DVf+ETsMheR8kroGGF+U5mjknVWXxtkOM+3E7X67czYn0TNrVq8p15zfg6ti6VAm3dZpN8TrrXkOliSUCE1DHDroG57DK8ONb8NUoCI+EFv1caaHZ5VAhosBTHE5J47PlOxn34y+s2/0rlSoEs+DhnkRWsrWZTfHxy5rFxhigUo6ZV2Kvhyp1YP2X8PP/YOV4t6jOn9e7nkeZma7KKZdqFUO5sVsjhndtyModh/kp4WB2Evj7F2tpUbsKQ+LPs4V0TImwRGDM2ahQCVr3d7eMdLcG8+4VJ7ufjv8dpCa7huZWv4HIUydMFBE6nhdJx/MiAUjPyGTZL0mEBAcxddlOHpm8kpQ018BtC+mY4mJVQ8YUp3nPwerJsG+de1ynI1zwB4i7vsCnpWdk0v35uXkupFM3MpwFD9tCOqboznZhGmPMmej+INy9EO5dCr3/BsFhkOytd5B6FL75K2xf6CbMyyEkOCjfhXR2JR1n5MTlfL16N8dSfZrqy5gCWYnAmJKm6sYobFsA7w+AzDSIiHbrKLTqD026Q0hY9tKa0RziXxVe5Z7U+9hHJBVDgwkLDSLpWBphIUFc0rwm/7wujogwq+k1BbMSgTGlRdZAtYYXwoNb4Np3oNHFsPpT+Hgw7P8ZgL9eWp2aoSe4L2QK58sG7g2ZTMXQYJ4e2J7Fj/Zi3G1d+V2XBqSkZVCpQjAAr83ZxMeLtgfqlZkyyr5CGBNI4VWh/SB3Sz/hGptj2gHQL/Ft+gSPQ8kkCGVoyDyi+/6VK7yG4m5No+jWNOqU0323cT91Iyvyuy4NUFXemLeFC5tG0aF+NcSH6TJM+WSJwJjSIiQMmvY8+Tj+FoJ2L4e9qwEII40rfroZui3L9xTjbu9KWobrZfTLwRRe+GYDGZlK7arh9G4TQ5+2MXRpHGXTZZtTWCIwprSqVg8ObDp1W9J2OLLXtSmMG+qql1oPgKim2YeEeiunNYiqxOJHezF7fSLfrN3DpCU7GLtwG1XCQ+jZqhZ92tSme8toKlvbQrlnjcXGlFZfjIRlY0+dHym4gltlrfuDMO462OWVDmLaQeurIPZ3p41VyHI8LYPvNu7nm7V7mLkukYNHU6kQHMSUuy+kbd1qqKpVH53jbGSxMWVNQZPkVakNt891JYR1X8C6z2HuM1A/3iWCpO1uZba6cdmN0+GhwfRqE0OvNjFkZCpLtx9i1rpEWsRUAeDZrzewbPshxt3W1dZQKGcsERhTWt3xXeHHRDaAbne525E9UMlrPF78Lnw3Bqqd50oKra+C87q4dZ1xayic36gG5zc6OV1Gw6hKnEjPyE4CIycuJ6ZqOH3axNCxfqQlh3OYVQ0Zcy46dhB+/hrWfg6bZ0PGCaje2A1sy2Puo9zSMzK5+b2f+GHzAdIzlVpVwujVJoY+bWLo1jSKsJDgEngRxt9s9lFjyqsTR2DjN66RuZu3ttR7/V1poc0AaNIj3/WaDx9LY86GRGas3cvcDYkcTc2gclgI3VtG06dNDD1a1aKqN322TZBX+lkiMMY46anw+b2w4X9w4jBUqAzN+8AFt0PDbvk+7XhaBj9sPsA3a/cwY+1e9ien8vLQWK6Jq8fHi7bxt2lrOZ5+cgW4rMFvlgxKD0sExphTpadCwnxYN81No93rCTcZ3q+7YctcaNkPKlbP86kZmcryXw7RIqYKVcJD6fjEdA6nnD7vUb3Iinz/cM88zmACwRKBMSZ/mRnuFlIBFv8XvrgfgkKg8aVunEKrK6FyrXyf3vjhL/Nd07Z7i2ha16lKm7pVaVOnKo1rRhBsDc8BYd1HjTH5CwrO7lFEpxFQuyOs+8w1Nn/xJ7cK2wObXAkhPdUljBzqRlbMc8rsiqHB7DtyggWbt5CW4VJFeGgQLWtXpVuTKB6+ohUAaRmZ2QPhTMkrsUQgIv8F+gOJqtrO21YDmAA0AhKAIap6qKRiMsbkISgI6nd2t15PwN41sHPJyWqicUPh+GGvW6ob1fxA35Y8MnkVldP2Z8+Umhwald1GkJqeyabEZNbu/pW1u35l3e5f+eXQsexLXvnKt8Q3qsE/ftsecHMmNatVmZiqYTbIrQSUWNWQiFwKJAMf5EgEzwEHVfUZEXkYqK6qDxV2LqsaMiaAFrzqFtvZtdQ9rtUWut7JVOlJ5hcjuSZ9OlNC+hLcf4xPDcWqyr/nbqZhVCX6d6jL4ZQ0Oj7xDQA1IirQuk4V2nhVS63rVKVpdGUrPZyhUtFGICKNgC9yJIINwGWqultE6gBzVbVlYeexRGBMKZD0C6z/wlUftewHHYbCPzu4WVSDK8A9S6B63tNdFCQ1PZMVO5JY55Ue1u7+lQ17jnDC65FUITiIFrUrc/dlzbiifR3SMjI5lppBtYqh/n6F55zS2kYQo6q7vft7gJhABmOMKYLI86Drne6mCl/++eRKaxmp8EosNL4EGnd3cyBVqe3TaSuEBJ026jk9I5Ot+4+6qiUvQWTNoLpyx2GufX0B799yAd1bRLP9wDHW7fmVNnWqUr96xXyrlmzcw0mBTgTZVFVFJN/iiYjcDtwO0KBB0b9lGGOKUfJeWP4RZObqQnp4F8x6Alpe4RJBwveQuBaaXAZRzU4u0lOIkOAgmsdUoXlMFa6OPfXDOqZqGA/2a0nrOm7OpG/W7uHvX7o1oquEh7geS1m3ulVpVqsyX6/ewyOTV5GS5hLXzqQUHpm8CqBcJoNAJ4K9IlInR9VQYn4HquqbwJvgqoZKKkBjjA/mPQeaeeq2oGBXIrj5SzdtNrjxCgtfc/er1HXLcjbu7qqVfJj6Ii/1q1firsuaZT++vktDOjesfkrD9ISffsn+0A8OEgRIzzz1YyQlLYPnp2+wRBAAnwM3Ac94Pz8LbDjGmDNS0EypOccf9H0Kzr8Vts6DLfPg5+mwfSHEDnP7l30I4ZFu6c6KkWcUSsUKwcQ1qE5cg5OD4TIylW0HjrJu9xHW7j7Ma3M25/ncXUkpvDl/MwePptGkZgRNoiNoXDOCGhEVzuneSyXZa2gccBlQE9gLPA5MBSYCDYBtuO6jBws7lzUWG3OOyMyE5D1Qta5rZ/hnBzeFtgS5KbQbd3eD2eqf1r55Vi56Znae4x7qRVakXb2qzF6fmD3uAaBqeAhNoivTpKZLDO3qVaNHq/wH2JVWAW8sVtVh+ey6vKRiMMaUMkFBLgmAay+4ZwnsXOymuNgyDxa84noh1Y+HjDT44V/Q6FKoG3tyANwZyBr3kFVdBG7wW1aDcXpGJjuTUtiy/yhb9h1l6/5ktu4/yg9bDjB52U66NYnKTgQ3vL2ITg0iGdnHdXj8ftN+6levSP3qlcrMCOpAVw0ZY8xJIRXc8psNL4Qef3Ezp6afcPv2roGZo939sGqu+qjJZW4GVR97JGXJagfIr9dQSHAQDaMiaBgVQY9cHdqPpabza455lepXr0jNKmHZ+65/exHgurk2iKpEY6+KyZUmKtMkOoKoUlbVZHMNGWPKjuRE2DrflRi2znPVSDd94Rql966B3StdA3RWKaOEpWVksvyXJLbuO8rm/cls3XeUrfuPsu3AMVIzTjam33d5c0b2bsHRE+m8891WftO+Ns1qVcn3vP7q6hrwqiFjjDlrlWtB+0HuBnBwK1T1PhDXTIX5z7n7Uc1dQmhyGbToB8ElM9gsNPj0MRDgGqt3Hkphy/5ktuw7SmyDSAC2HTjGmBk/0yKmCs1qVWHB5v38eeIKGnttEY1rRpD463He/2Fb9oC64ujqaiUCY8y5ITMTEte4toUtc2HbAggOgQe3uvaE9V9BaEVo0NX9BLe856SbYdB7UCUw41lTUjMQcWtKr955mP9+v5Ut+46yZV8yvx4/fWrvLGcyxXepmGLCXywRGGMKlZ4Kh7ZCtFfJ/+9ubjBbcBicd4ErMexZ5dZj6Hwz9B8T2HhzUVUOHUuj85Mz8pziW4Ctz1xZpHNa1ZAxpnwJqXAyCQDc+g1s+8EbwzAXZv/ddVPVTDcqutp50KIPRLc+48Ft/iQi1IiokO8U33UjK/rvWlYiMMaUS1PuhFWfQGYaBIW6nwCVanpzJF0KLa6AqnUCGubUZTvz7Op6JsuA5lciCHzaM8aYknZkD6yZfPLDPzMNQsKg79PQrJcb7fzF/fDLQrf/UAIsHweHd5Z4qNfE1ePpge2pF1kRwbUN+HstaKsaMsaUP3nNjaQKBzbBwP949zefbEDe8DV87S2VUqOpKy006e56JIX6r4omP9fE1SvWOZAsERhjyp+C5kYCN8q55smJ7LjgdjeAbes8N45h1SRYNhYe2ub2b5rpRj43vAjCq5bMa/AjSwTGmPLnju+KdnxQENRu527d7oaMdNj/M4RVdvu/f8UlCQn25ki6FJpd7pJHGWCJwBhjiio4BGLanHz8u4mw4ydXWtg6382RtHv5yUTw0zsQ0xbqdnK9mUoZSwTGGHO2QsO9nkaXAI/CiWQ4tt/tO/4rfPUAaAaERkDDbq7E0Ko/RDUNaNhZLBEYY4y/hVU+WW0UXhUe2AQJ350sMcx4DEIruURwZC+s/cwlh+iWPq/a5k+WCIwxprhVquFmSW0zwD0+ssd1VwXY9j387wF3v3KMSwiNL4U2V0N4tZPnKMbpMCwRGGNMScs5bXa7gVCv08nSwtb5bqBb08tdItg63yWBzbPd+IZ5z/p9OgxLBMYYE2jVG7lbpxtPjmGo5o0bWPYhrJxw8tjlH0H3h/xaKrCRxcYYU5rkHsNwzevQ+mrXNRXcQLh5z/r1kpYIjDGmNDu6DzZOd72OwA18W/6Ra2T2E0sExhhTmuU5HYZ/SwWWCIwxpjQrbDoMP7DGYmOMKc2KOh3GGbASgTHGlHOWCIwxppyzRGCMMeWcJQJjjCnnLBEYY0w5VyYXrxeRfcC2M3x6TWC/H8MpbmUpXou1+JSleMtSrFC24j3bWBuqanTujWUyEZwNEVmsqvGBjsNXZSlei7X4lKV4y1KsULbiLa5YrWrIGGPKOUsExhhTzpXHRPBmoAMoorIUr8VafMpSvGUpVihb8RZLrOWujcAYY8ypymOJwBhjTA6WCIwxppwrN4lARP4rIokisjrQsRRGRM4TkTkislZE1ojIHwMdU0FEJFxEfhSRFV68TwQ6psKISLCILBORLwIdS2FEJEFEVonIchFZHOh4CiIikSIySUTWi8g6EekW6JjyIyItvfc06/ariPwp0HHlR0Tu9/6/VovIOBEJ99u5y0sbgYhcCiQDH6hqu0DHUxARqQPUUdWlIlIFWAJco6prAxxankREgAhVTRaRUOA74I+qujDAoeVLREYC8UBVVe0f6HgKIiIJQLyqlvpBTyLyPvCtqr4tIhWASqqaFOCwCiUiwcBOoIuqnulg1WIjIvVw/1dtVDVFRCYCX6nqe/44f7kpEajqfOBgoOPwharuVtWl3v0jwDqgXmCjyp86yd7DUO9War9hiEh94Erg7UDHci4RkWrApcA7AKqaWhaSgOdyYHNpTAI5hAAVRSQEqATs8teJy00iKKtEpBEQBywKcCgF8qpalgOJwAxVLc3xvgw8CGQWclxpocA3IrJERG4PdDAFaAzsA971qt3eFpGIQAflo+uAcYEOIj+quhN4AdgO7AYOq+o3/jq/JYJSTEQqA58Cf1LVXwMdT0FUNUNVY4H6wAUiUiqr30SkP5CoqksCHUsRXKyqnYArgLu9as7SKAToBLyuqnHAUeDhwIZUOK8KawDwSaBjyY+IVAeuxiXbukCEiNzgr/NbIiilvLr2T4GPVHVyoOPxlVcVMAfoF+BQ8nMRMMCrdx8P9BSRDwMbUsG8b4OoaiIwBbggsBHlawewI0dpcBIuMZR2VwBLVXVvoAMpQC9gq6ruU9U0YDJwob9ObomgFPIaX98B1qnqmEDHUxgRiRaRSO9+RaA3sD6gQeVDVR9R1fqq2ghXHTBbVf32zcrfRCTC6zCAV83SByiVPd9UdQ/wi4i09DZdDpTKDg65DKMUVwt5tgNdRaSS9/lwOa7t0C/KTSIQkXHAD0BLEdkhIrcGOqYCXAQMx31bzera9ptAB1WAOsAcEVkJ/IRrIyj13TLLiBjgOxFZAfwIfKmqXwc4poLcC3zk/S3EAv8IbDgF85Jrb9w37FLLK2VNApYCq3Cf3X6bbqLcdB81xhiTt3JTIjDGGJM3SwTGGFPOWSIwxphyzhKBMcaUc5YIjDGmnLNEUI6JyFwR+VeArq0iMigQ1/ZFaY/PH7xZLEf7cNwcEbmxBELK69oF/h68WW+vLcmYzkWWCM5R3iCvf3tTGJ8Qkb0iMktEeuc4bCDwSKBi9DcR6eR9cFySz/4JIrKgpOMqSH4fdCLyXmmYIltErgTOAz7KsS3Bi1tFJMWbcvoBb6BTSXsSeEZE7LPsLNibd+76FDcVwa1AC6A/8D8gKusAVT3ozW5a5ohISO4PHm/G1uXALXkcHwVcg804WlR/BN5T1Yxc2/+GG0jYGjcZ2j+AQEyI9xVQBTdNhDlDlgjOQd50D5cAD6vqLFXdpqo/qeoLqjo+x3GnVA153/T+T0T+4y3SsUNEHsh17hYiMk9EjovIBhH5jYgki8gIb38j75tifK7nFVbEf8Y7X4oXx3M5F94QkdFeVcYIEdkMnADymtnybWCwN2FfTjd4z5kgIv1E5FsROSQiB0Vkuoi0LiA2n16TiNQTkfHeeQ+JyJci0jy/8xaFiAwUkZXe+3PQ+x3E5Nh/lbjZSY+LyFYReUrcZGpZ+2uJyGfe87eJyGnJMo9rRuPmuJmWx+4jqrpHVRNU9W1gJW76i6znNvWut0dEjorIUnET/uU8f6F/b3nE9JCI7BeRruAmO8Qlg2GFvR6TP0sE56Zk7zZAir6K0f24IeydgGeB58RbZcorfk8B0oGuwAjgcSDMDzEfxX2Tbw3chZsH6NFcxzQGfgcMBjoCx/M4z0dAMDA01/ZbgQmqehSXQF7GlZguAw4D03J+cBaViFTCTbZ3HOgOdMNNFzzT23fGRKQ2boK893Hvz6XA2Bz7++Je97+Atrj3cRCnTu/wHtAM98F+DXAj0KiQS1+MS575zm0kzmVeXGk5dlXGlUB7435XnwKTRaRVrlPk+/eWx3VewE1h0T3Xokc/4t5zc6ZU1W7n4A24FrcQz3HcHEsv4FZfynnMXOBfOR4nAONyHbMR+D/vfl9cEqiXY/+FuPnyR3iPG3mP43OdR4FB+T3OI/47gE05Ho/GfdDE+PDaPwQW5Hh8vne9LvkcHwFk4KZ7Pi0+X14T7sN3I960Ld62YOAAMKSAWPN8H3Af3F949zt5xzXM5xzzgb/m2nYN7suA4KoGFbgox/6G3mseXUBsfwK25bE9AZcgkoFU79wpwIWF/F4WZv0t+fL3luP9GQq8C/yc13uAm0I6Ewgp6f+zc+VmJYJzlKp+ipu3/CrcN7MLgYUi8pdCnroy1+NdQC3vfitgl3rTInt+wg8LvIjIIBH5zqtKSAZeAhrkOmyH+jZV8NtAtxzfPm8BVqs3PbJXbfGxiGwWkV+BvbjSce7rFUVnXInliFdVlowraVQHmp7FeQFWADOB1SLyqYjc6VXb5Lz2o1nX9a79MS7B1cZ9W8/EfXMGQN1KXIWtcFWRvEtdAGNwk8p1x5WEnlDV7IZ4cbOmPidu3e1DXkzxnP4eF/T3luUFXMntYs17BbEUXMLz2xq+5Y0lgnOYqh5X1Rmq+jdVvRA3tfXoQqpA0nI9Vor2d5KVFLIbcsWtrZAvr753PDAdl7jigP/DLXmZ01EfY5gHbAJuETct9jC85RM9XwDRwB+ALt710oH83hdfXlMQrqE6NtetBfCfAmI9AlTLY3skLpGgrh68j3dbiavm2igiHXNc+4lc1+0ANMetGJalqDNM7sclsrwcUNVNqvoDrvQ5SkR65Nj/Aq4K76+4ZBGLS0S532Nf/t5m4BJafjPw1gCO68nlUk0RhQQ6AFOi1uJ+5+G4In1RrQfqikhdVc36NhnPqf+4WR88dXJsiy3kvBcBO1X1yawNItLwDOID3BrKIvJfXI+X9bhvtmO980bhSjZ3qeocb1snCv5f8OU1LcUlnP1atHV6N+C+0WcnKnELqXfEVYdkvyZcFd8PIvI3YA2uymSFd+1WqroprwuIyHrc7+gCYIG3rQGuxFiQZUC0iNRU1f35HaSqh8R1OnhJROK8WC8GPvBKpnhtVU1x1TtF9RVumuhPRERV9f1c+9vh3gNzhqxEcA4SkSgRmS0iN4hIBxFpLCKDcev0ztIzX/ZyBu6D630R6eh9kx+D+zbtKnRVU3B1wQ+JSFsRuRD37bAgPwP1ROR6EWkiIndy9r1A3gdqeteeqqoHvO2HcN90bxORZiLSHXjDew158vE1fYSrYvpMRLp77/mlIvJiIT2HxuBKLneL65EVi5tnvob3ExHp6vWuOd/7AB+A69uftejL34DficjfRKSdiLTyqtqe8+LfAHwN/EdEunnXeA9XpVKQZbg1qC8u5DiAfwMtcaUAcL/T34ob29Ee125zxlU36ta3GAy8IacPbrsE9/rMGbJEcG5Kxn1w/RFXTbIG14PkY07vTeMzVc0EfovrJfQj7sP2KVwSyFmXnNU18Sdctcj/FXLeacDzuJ48K3E9TR470zi9c+7CfZOsTo6xA95rGIqrOlkNvIarvjhRyCkLfE2qegzXm2cLbu3b9bj3pzou+eQX5zjgZu+2GPeBVhu4RN2KX+CqiC7CVWltBF4EnlTVD71zTAeuBHrgfi8/4tYK3p7jUiOArcBsXHfQj3GNtfnyqqT+C1xf0HHesYm4Utdor3fZSFwS+RbXRrXQu3/GvGQwBJfQbgTXZRfX/vVuQc81BbOFacxZ8eqpl+N61JSlBeGND0SkFq7kcb6qbg10PLmJyPNANVUNxGC2c4a1EZgiEZHf4hptN+K6VY7hZD21OceoaqK4wWcNcCWK0iaRwqseTSGsRGCKxCuS/x+ujvoQbizC/T526zTGlEKWCIwxppyzxmJjjCnnLBEYY0w5Z4nAGGPKOUsExhhTzlkiMMaYcu7/AXkI/NK2++0sAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "singular_value = singular_value_list[-1]\n",
    "err_subfull, err_local, err_SVD = [], [], []\n",
    "U, D, V_dagger = np.linalg.svd(M, full_matrices=True)\n",
    "\n",
    "# 计算 Frobenius-norm 误差\n",
    "for i in range(T):\n",
    "    lowrank_mat = np.matrix(U[:, :i]) * np.diag(D[:i])* np.matrix(V_dagger[:i, :])\n",
    "    recons_mat = np.matrix(U_learned[:, :i]) * np.diag(singular_value[:i])* np.matrix(V_dagger_learned[:i, :])\n",
    "    err_local.append(norm(lowrank_mat - recons_mat)) \n",
    "    err_subfull.append(norm(M_err - recons_mat))\n",
    "    err_SVD.append(norm(M_err- lowrank_mat))\n",
    "\n",
    "\n",
    "# 画图 \n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(list(range(1, T+1)), err_subfull, \"o-.\", \n",
    "        label = 'Reconstruction via VQSVD')\n",
    "ax.plot(list(range(1, T+1)), err_SVD, \"^--\", \n",
    "        label='Reconstruction via SVD')\n",
    "plt.xlabel('Singular Value Used (Rank)', fontsize = 14)\n",
    "plt.ylabel('Norm Distance', fontsize = 14)\n",
    "leg = plt.legend(frameon=True)\n",
    "leg.get_frame().set_edgecolor('k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### 案例2：图像压缩\n",
    "\n",
    "为了做图像处理，我们先引入必要的 package。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:17.225100Z",
     "start_time": "2021-03-09T03:50:17.151352Z"
    }
   },
   "outputs": [],
   "source": [
    "# 图像处理包 PIL\n",
    "from PIL import Image\n",
    "\n",
    "# 打开提前准备好的图片\n",
    "img = Image.open('./figures/MNIST_32.png')\n",
    "imgmat = np.array(list(img.getdata(band=0)), float)\n",
    "imgmat.shape = (img.size[1], img.size[0])\n",
    "imgmat = np.matrix(imgmat)/255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:18.337792Z",
     "start_time": "2021-03-09T03:50:17.231211Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUsUlEQVR4nO3dbWxcVXoH8P8fx46NHSW2A4mdOC+khGAqNomsQEVAdOkugS8EtOLlA2Il1KwqkIq0+wFRbZdWVctWBUSliioUtNmKQtgFRFqhdilaKUVasWvSxCREJCE4L7bjJMRJnMRxYufph7nZOtF9jsczd2Zsn/9Psjw+z5yZx9fz+I7v8TmHZgYRmf6uqXQCIlIeKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFil6KR/D7JUZJnxnzcXem85EozKp2ATBu/MbO1lU5CfDqzT3Mku0n+iGQXyVMkN5OsrXReUn4q9jg8DGAdgKUAbgXw/bQ7kVxL8mTgI3TmXkXyOMk9JH9MUu8aJxn9QOLwj2bWCwAk/x3AyrQ7mdknAOYU8PhbAfwhgAMAbgGwGcAIgL8r4LGkRHRmj8ORMbfPAWjI8sHNbL+ZfW1ml8zscwB/DeB7WT6HFE/FLr9H8s6rrqhf/XFnng9lAFjKXGXi9DZefs/M/gcFnPVJ3gdgm5n1k1wB4McAfpF1flIcndklC/cA6CJ5FsCHAN4D8LeVTUmuRi1eIRIHndlFIqFiF4mEil0kEip2kUiUdeitsbHRWltby/mUmbrmmvTfjV47AJD+cHPo4ujo6GhBsUIuuBaafyh26dKlCbUD4dxnzPBfqlVVVW7MU2gek11vby8GBgZSfzBFFTvJdQBeAVAF4F/M7IXQ/VtbW7F58+ZinrKiZs6cmdpeX1/v9qmpqXFjFy5ccGOnTp1yY6dPn3Zjw8PDqe2hgq6rq3Nj3vcMhAtwaGgotX1wcNDtMzIy4saam5sLinmFe+7cObfPxYsX3dhk98gjj7ixgt/Gk6wC8E8A7gPQDuAxku2FPp6IlFYxf7OvAbAv+b/oCwDeBvBANmmJSNaKKfYFAA6N+fpw0nYFkhtIdpLsHBgYKOLpRKQYJb8ab2YbzazDzDoaGxtL/XQi4iim2HsAtI35emHSJiKTUDHF/jsAN5JcSrIGwKMAtmSTlohkreChNzMbIfk0gP9CbujtDTPblVlmIpKposbZzexD5KY0isgkp3+XFYmEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4mEil0kEip2kUio2EUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEUTvCkOwGMAhgFMCImXVkkZSIZK+oYk/8sZkdz+BxRKSE9DZeJBLFFrsB+BXJz0huSLsDyQ0kO0l2DgwMFPl0IlKoYot9rZmtBnAfgKdI3nX1Hcxso5l1mFlHY2NjkU8nIoUqqtjNrCf5fBTA+wDWZJGUiGSv4GInWU9y1uXbAL4LYGdWiYlItoq5Gj8PwPskLz/Ov5nZf2aSlYhkruBiN7P9AL6VYS4iUkIaehOJhIpdJBIqdpFIqNhFIpHF/8ZH49SpU6ntPT09bp/z58+7sZkzZ7qx2bNnu7Frr73WjY2Ojqa29/f3u33OnDkz4ccDgAsXLrixwcHBCedx+PBhN1ZVVeXG2tvb3diKFStS2xcuXOj2qaurc2NTmc7sIpFQsYtEQsUuEgkVu0gkVOwikSj71XgzK/dTZubIkSOp7bt27XL7dHd3u7Ha2lo31tTUVFA/78p6b2+v2+fYsWNuzBuBAIChoSE35uWYzKVIdfDgQTe2f/9+N7ZkyRI3tn79+tT2e++91+2zYMECNzaV6cwuEgkVu0gkVOwikVCxi0RCxS4SCRW7SCQ0EWYCzp07l9p+6NAht09XV5cbO336tBsLTQoJDYd5Q17XXXed26ehocGNhSbJjIyMuLGlS5dOqH28PELHOBTzJt4MDw+7faYrndlFIqFiF4mEil0kEip2kUio2EUioWIXiURZh95IoqamppxPmam2trbU9jvvvNPts2jRIjfmzaIDgJ07/Z20+vr63Nj8+fNT2++44w63z2233TbhxwOAa67xzxXe8OCePXvcPlu3bnVjt9xyixsLDeetXbs2tb2lpcXtM5Vfo6FZheOe2Um+QfIoyZ1j2ppIfkRyb/JZ27OKTHL5vI3/GYB1V7U9C+BjM7sRwMfJ1yIyiY1b7Ga2FcCJq5ofALApub0JwPps0xKRrBV6gW6emV3+w/EIcju6piK5gWQnyc6BgYECn05EilX01XjLrTPlrjVlZhvNrMPMOhob9ae9SKUUWuz9JFsAIPl8NLuURKQUCh162wLgCQAvJJ8/yLfjVF5w0tt2KbRA4Zw5c9zY6tWr3dhDDz2Ud15jedtNXbx40e0TynHu3LlurL6+3o1dunQptf3Eiasv//y/vXv3urHQDLvQ8Oa8eel/YVZXV7t9pvJrNCSfobe3APwGwE0kD5N8Erki/w7JvQD+JPlaRCaxcc/sZvaYE7on41xEpIT077IikVCxi0RCxS4SCRW7SCS019sEeDOKQsM4M2fOdGOhPdtCC0SGZmWdPXs2tT20SGUoxxkz/JfI6OioG/P2gQstshmKeUN5QHj2XSjmmcqv0RCd2UUioWIXiYSKXSQSKnaRSKjYRSKhYheJhIbeMhAa3gkNXVVVVbmx0FBTSF1dXWp7KMdQHqHhtdBiJL29vanthexTB/gzDoHwrL3QsKhnOr5GAZ3ZRaKhYheJhIpdJBIqdpFIqNhFIlHWq/FmVvBV5sksdDU7NGkl1C+05lroGHpX42fNmuX2CV1xD109P3DggBvbt2/fhB8vtN5daJ0/b505wD8eIVP5NRoaSdCZXSQSKnaRSKjYRSKhYheJhIpdJBIqdpFIaCLMBHi5e2vTAeHhtVC/0HZNoWE5b+JHQ0OD28fbMgoADh8+7Ma++uorN7Z///7U9tA6c83NzW5syZIlBfULra/nmcqv0ZB8tn96g+RRkjvHtD1Psofk9uTj/tKmKSLFyudt/M8ArEtpf9nMViYfH2ablohkbdxiN7OtAPytN0VkSijmAt3TJLuSt/mN3p1IbiDZSbIztNiBiJRWocX+KoBlAFYC6APwondHM9toZh1m1tHY6P5OEJESK6jYzazfzEbN7BKA1wCsyTYtEclaQUNvJFvMrC/58kEAO0P3ny68obLQ8FooFhriCQ2HhdaT84blQn2Gh4fdWE9Pjxv74osv3Jg39BZaS2758uVubNWqVW4sNFvOm3VYyLZQU924xU7yLQB3A5hL8jCAnwC4m+RKAAagG8APSpeiiGRh3GI3s8dSml8vQS4iUkLxvZcRiZSKXSQSKnaRSKjYRSJR9llvoZlek503VBYaQgt9v6F+oZlthSxGGVpEcWhoyI0dP37cjYWG5fr7+1Pb29ra3D5NTU1u7Prrr3djoUUlve879HOZyq/REJ3ZRSKhYheJhIpdJBIqdpFIqNhFIqFiF4lEWYfeSE7p2UaFDL2FhI5FocfJG0Y7duyY2+fQoUNurLu7240dOXLEjV24cCG1PbTnXGtrqxsL7ecWGooMDSt6pvJrNDRsOHW/KxGZEBW7SCRU7CKRULGLRELFLhIJTYTJQKETYUJXfUPbFs2Y4f/Yzpw5M6F2APjyyy/dWOhqfGgrp9ra2tT20ISW+fPnu7HZs2cXlEdoLT/PdHyNAjqzi0RDxS4SCRW7SCRU7CKRULGLRELFLhKJfHaEaQPwcwDzkNsBZqOZvUKyCcBmAEuQ2xXmYTOLcpvW0NBbKBbaGsrbtggIDw2dPHkytf3rr792++zZs8eNhSa7hIYAvSG2lpYWt09okkzI6OioGytkDbrpKp8z+wiAH5pZO4DbATxFsh3AswA+NrMbAXycfC0ik9S4xW5mfWa2Lbk9CGA3gAUAHgCwKbnbJgDrS5SjiGRgQn+zk1wCYBWATwHMG7OT6xHk3uaLyCSVd7GTbADwLoBnzOyK/0+03B+mqX+cktxAspNk54kTJ4pKVkQKl1exk6xGrtDfNLP3kuZ+ki1JvAXA0bS+ZrbRzDrMrCO0CYCIlNa4xc7cZcvXAew2s5fGhLYAeCK5/QSAD7JPT0Syks+stzsAPA7gc5Lbk7bnALwA4B2STwI4AODhkmQ4iXjDaKGtlUJCQ1ehWV7e+m6AP0utq6vL7bNt2zY3dvRo6hs2AOHtmtrb21Pbb7rpJrdPQ0ODGwvNXhseHnZj3rBcaNhzuhq32M3sEwDeoOQ92aYjIqWi/6ATiYSKXSQSKnaRSKjYRSKhYheJRNkXnCx0q6TJwBtiCw29hb7fUL/Q0NDFixfdWE9PT2r7jh073D6hBSdDM9FWrFjhxtauXZvafuutt7p9QotserP5AGBwcNCNeUKLfU7l12iIzuwikVCxi0RCxS4SCRW7SCRU7CKRULGLRKKsQ29mVvAMscnAW6QwtHhhKBYa4hkZGXFjQ0NDbswbhgoNT4XyaG5udmPLly93Y97stnnz/AWNQnu2hfaqC80CrK6uTm0vdEh0sgt9Xzqzi0RCxS4SCRW7SCRU7CKRULGLRKLsE2Gmo0InVYSuIn/zzTdurK+vz42dOnUqtb22ttbt09ra6sZuuOEGN7Zs2bIJP2Z9fb3bx8sdCE/+CfF+Ntr+SUSmLRW7SCRU7CKRULGLRELFLhIJFbtIJMYdeiPZBuDnyG3JbAA2mtkrJJ8H8KcAjiV3fc7MPixVopNBIRNhvO2HgPDEj+PHj7sxb4snAOjv709tb2xsdPssXrzYjd18881urK2tzY15a9eF1tYrdEJRiDf0Fhouna7yGWcfAfBDM9tGchaAz0h+lMReNrN/KF16IpKVfPZ66wPQl9weJLkbwIJSJyYi2ZrQexmSSwCsAvBp0vQ0yS6Sb5D03yeKSMXlXewkGwC8C+AZMzsN4FUAywCsRO7M/6LTbwPJTpKdAwMDxWcsIgXJq9hJViNX6G+a2XsAYGb9ZjZqZpcAvAZgTVpfM9toZh1m1hG6SCQipTVusTN3GfR1ALvN7KUx7S1j7vYggJ3ZpyciWcnnavwdAB4H8DnJ7UnbcwAeI7kSueG4bgA/KEF+k4o3XFPoemahLY12797txvbu3evGvJl0CxcudPuE1pJbtGiRG2toaHBj58+fT20PDUWGhsNCs/ZCM+JmzEh/icc46y2fq/GfAEg7MtN6TF1kuonvPwtEIqViF4mEil0kEip2kUio2EUioQUnSyw0LOcNTwHAoUOH3NjBgwfdWFNTU2p7aIZae3u7G5szZ44bC/H+W9IbChtPaOgtNIzmzbILDfNN5e2fQnRmF4mEil0kEip2kUio2EUioWIXiYSKXSQSGnrLQGgmV2g/t8HBQTfW29vrxg4cOODGvGGo6upqt8/s2bPdWF1dnRsbGRlxY96wYmgILRQLLVQZGirzYjEuOBnfdywSKRW7SCRU7CKRULGLRELFLhIJFbtIJDT0NgHeME5o6G14eNiNnT171o0dO3asoNjcuXNT20OzzUILR4b6hb43bxHI0JBXTU2NGwsNr4UWnPSGIjX0JiLTlopdJBIqdpFIqNhFIqFiF4nEuFfjSdYC2ApgZnL/X5rZT0guBfA2gGYAnwF43Mz8WR+RKsU2Q6FJId5kktAV9/r6ejcWWkMvNAoxNDTkxjyFXo2PcSunQuRzZh8G8G0z+xZy2zOvI3k7gJ8CeNnM/gDAAIAnS5aliBRt3GK3nDPJl9XJhwH4NoBfJu2bAKwvRYIiko1892evSnZwPQrgIwBfAThpZpcnNB8GsKAkGYpIJvIqdjMbNbOVABYCWANgRb5PQHIDyU6Snd5a4iJSehO6Gm9mJwH8GsAfAZhD8vIFvoUAepw+G82sw8w6Ghsbi8lVRIowbrGTvI7knOR2HYDvANiNXNF/L7nbEwA+KFGOIpKBfCbCtADYRLIKuV8O75jZf5D8AsDbJP8GwP8CeL2EeU4K3uSJ0NppoeGp5uZmN7Z48eKCHnPRokWp7aFtnELr04WeKzTk5Q2VhSbPhGKhCTmhdfK8YcpQ7qHveSobt9jNrAvAqpT2/cj9/S4iU4D+g04kEip2kUio2EUioWIXiYSKXSQSDM1qyvzJyGMALu9dNBfA8bI9uU95XEl5XGmq5bHYzK5LC5S12K94YrLTzDoq8uTKQ3lEmIfexotEQsUuEolKFvvGCj73WMrjSsrjStMmj4r9zS4i5aW38SKRULGLRKIixU5yHckvSe4j+Wwlckjy6Cb5OcntJDvL+LxvkDxKcueYtiaSH5Hcm3wu+UofTh7Pk+xJjsl2kveXIY82kr8m+QXJXST/PGkv6zEJ5FHWY0KyluRvSe5I8virpH0pyU+TutlM0l+ON42ZlfUDQBVya9jdAKAGwA4A7eXOI8mlG8DcCjzvXQBWA9g5pu3vATyb3H4WwE8rlMfzAH5U5uPRAmB1cnsWgD0A2st9TAJ5lPWYACCAhuR2NYBPAdwO4B0Ajybt/wzgzybyuJU4s68BsM/M9ltunfm3ATxQgTwqxsy2AjhxVfMDyK3SC5RptV4nj7Izsz4z25bcHkRuJaQFKPMxCeRRVpaT+YrOlSj2BQAOjfm6kivTGoBfkfyM5IYK5XDZPDPrS24fATCvgrk8TbIreZtf1oUDSS5BbrGUT1HBY3JVHkCZj0kpVnSO/QLdWjNbDeA+AE+RvKvSCQG53+zI/SKqhFcBLENuQ5A+AC+W64lJNgB4F8AzZnZ6bKycxyQlj7IfEytiRWdPJYq9B0DbmK/dlWlLzcx6ks9HAbyPyi6z1U+yBQCSz0crkYSZ9ScvtEsAXkOZjgnJauQK7E0zey9pLvsxScujUsckee6TmOCKzp5KFPvvANyYXFmsAfAogC3lToJkPclZl28D+C6AneFeJbUFuVV6gQqu1nu5uBIPogzHhLnVH18HsNvMXhoTKusx8fIo9zEp2YrO5brCeNXVxvuRu9L5FYC/qFAONyA3ErADwK5y5gHgLeTeDl5E7m+vJ5HbIPNjAHsB/DeApgrl8a8APgfQhVyxtZQhj7XIvUXvArA9+bi/3MckkEdZjwmAW5FbsbkLuV8sfznmNftbAPsA/ALAzIk8rv5dViQSsV+gE4mGil0kEip2kUio2EUioWIXiYSKXSQSKnaRSPwffyIxG32ku9UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUkklEQVR4nO3df4xV5Z3H8fe3yK/KYBGmgIgDAroCRSATUiq2rm4t2jTaZNNWN103caXp1mSbdP8w3WTrbvaPdrNt081u3NDV1DZdqbaauhvaqtQGTS064AhYrfwoWJAfQwUB5YcD3/3jHrID3u8zM+f+muH5vJLJ3Hm+99zzzJn7nXPv+d7neczdEZHz3/ta3QERaQ4lu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULLLoJnZfDP7hZkdMLP3fFDDzC42s8fM7G0z22lmt7ein3I2JbuU8S7wMHBnEP8P4CQwGfgL4D4zm9ekvknA9Am684uZ7QD+HfhLoAP4OXCHux9vwL5mA1vc3fq0XQgcBOa7+2tF2w+A3e5+T737IAOnM/v56TPAcmAmsAD4q2p3MrNlZnYo8bWsxL6vAHrPJHrhJUBn9ha7oNUdkIb4N3d/A8DM/gdYWO1O7v4s8IE673sccPictreAtjrvRwZJZ/bz094+t9+hkoDNchQYf07beOBIE/sgVSjZM2Zm15rZ0cTXtSUe9jXgAjOb06ftauDl+vRaytLL+Iy5+zOUOOubmQGjgVHFz2MqD+cn3P1tM3sU+Ccz+2sqbyFuAT5St45LKTqzSxkdwDH+/2x9DPhdn/jfAGOB/cBDwBfdXWf2FlPpTSQTOrOLZELJLpIJJbtIJpTsIploault0qRJ3tHR0cxdimRl586dHDhwwKrFakp2M1sOfAcYAfyXu389df+Ojg5+/etf17LLIalSdq4uVe1Ixco+ZpnHSynbx2Y9Xn+PWe99DXUf+Uj8cYbSL+PNbASVoYw3AXOB28xsbtnHE5HGquU9+xJgq7tvd/eTwCoqn5QSkSGolmSfBvyhz8+7irazmNkKM+sys66enp4adicitWj41Xh3X+nune7e2d7e3ujdiUiglmTfDUzv8/OlRZuIDEG1JPsLwBwzm2lmo4DPAY/Xp1siUm+lS2/u3mtmdwO/oFJ6e0Ajm0SGrprq7O6+Glhdp76ISAPp47IimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimVCyi2RCyS6SCSW7SCaU7CKZULKLZELJLpIJJbtIJpTsIplQsotkQskukgklu0gmlOwimahpRRgz2wEcAU4Bve7eWY9OiUj91ZTshT919wN1eBwRaSC9jBfJRK3J7sATZrbezFZUu4OZrTCzLjPr6unpqXF3IlJWrcm+zN0XAzcBXzKzj557B3df6e6d7t7Z3t5e4+5EpKyakt3ddxff9wOPAUvq0SkRqb/SyW5mF5pZ25nbwI3A5np1TETqq5ar8ZOBx8zszOP8t7v/vC69EpG6K53s7r4duLqOfRGRBlLpTSQTSnaRTCjZRTKhZBfJRD0+G5+Nt956a1DtAO97X/z/dObMmWHs+PHjYWzEiBFh7MiRI1Xbd+3aFW5z+vTpMHby5Mkwtn379jB29OjRqu2p45E6jtu2bQtjF1wQP43nz59ftX3ZsmXhNrNmzQpjw5nO7CKZULKLZELJLpIJJbtIJpTsIplo+tV4d2/2Lutmx44dVdufe+65cJt9+/aFsblz55bqx6lTpwa9v9///vfhNqm/yeHDh8PY66+/HsaKMRPvceGFF4bbpK7Gd3d3h7He3t4wduutt1Ztv+KKK8JtLr/88jA2nOnMLpIJJbtIJpTsIplQsotkQskukgklu0gmNBBmEKLS0Pr168Ntfvazn4Wxt99+O4yNGTMmjKVKb9GgkPHjx5faV2pAzrhx48LYpZdeWrV9ypQp4TYTJkwIY1HZE9Ilu2hG44kTJ4bbnK90ZhfJhJJdJBNKdpFMKNlFMqFkF8mEkl0kE00vvaXmIBvqFi9eXLU9NZIrNbpq3bp1YSw12iy1v6iP119/fbjN7Nmzw1iqZJdalfeiiy6q2j569Ohwm9WrV4exX/3qV2Fs3rx5YeyGG26o2p6a/284P0dT+v2tzOwBM9tvZpv7tF1sZk+a2Zbie1wgFZEhYSD/wr4HLD+n7R5gjbvPAdYUP4vIENZvsrv7WuDNc5pvAR4sbj8I3FrfbolIvZV9czLZ3fcUt/dSWdG1KjNbYWZdZtZ14MCBkrsTkVrVfCXCK3MahfMauftKd+90985JkybVujsRKalssu8zs6kAxff99euSiDRC2dLb48AdwNeL7z8d6IbDecLJqJy0aNGicJtUiee2224LY2XLP1FpKzVCLVUOS42wmz59ehiLlqhKTRy5atWqMBYtJwWwdOnSMHbZZZeFschwfo6mDKT09hDwHHClme0yszupJPnHzWwL8GfFzyIyhPV7Znf36PRT/dMKIjIknZ8fFRKR91Cyi2RCyS6SCSW7SCa01tsgROWw1ISNqUkU29rawtjYsWPD2MmTJ8NYdHxHjhwZbpNy4sSJMJYafRdNpvnqq6+G22zYsCGMRWvHQbq8+cEPfrBq+6hRo8JtTp8+HcaGM53ZRTKhZBfJhJJdJBNKdpFMKNlFMqFkF8mE1nobhN7e3qrtqbJQKhaNDIPyJcrU/iKpUlNq9F2qj9u2bava/vzzz4fbHDp0KIwtXLgwjM2dOzeMRaXP1O+l0puIDGtKdpFMKNlFMqFkF8mEkl0kE02/Gl/mavFQEV2NT129TQ1AueCC+PC/++67YSx1DKOrzKk+puaZSw3IiQa7AKxfv75q+7PPPhtukzpWN910Uxj70Ic+FMai/qeOx3B+jqbozC6SCSW7SCaU7CKZULKLZELJLpIJJbtIJjQQZhCiElWqjJMqr6UGY6TKYan506LHTJXyUlKDdfbs2RPGotLb1q1bw21SSzUtWbIkjLW3t4exSFRGhfTvPJwNZPmnB8xsv5lt7tN2r5ntNrPu4uvmxnZTRGo1kJfx3wOWV2n/trsvLL5W17dbIlJv/Sa7u68F3mxCX0SkgWq5QHe3mW0sXuaHk6Ob2Qoz6zKzrp6enhp2JyK1KJvs9wGzgIXAHuCb0R3dfaW7d7p7Z5kLKSJSH6WS3d33ufspdz8NfBeIL5WKyJBQqvRmZlPd/Uzd5dPA5tT9zxdlllBKlddS5bDUdqlRWVEZMFXKSy1fdezYsTD29NNPh7ForrmJEyeG23z2s58NYwsWLAhjKWXmDTxf9ZvsZvYQcB0wycx2AV8DrjOzhYADO4AvNK6LIlIP/Sa7u99Wpfn+BvRFRBpIH5cVyYSSXSQTSnaRTCjZRTKhUW+DEI2GKjt6LTXyKjVaLlU2ivZXtpT3xz/+MYx1d3eHsR07dlRtv+SSS8JtUss4tbW1hbHUcYyWqEodj/NVfr+xSKaU7CKZULKLZELJLpIJJbtIJpTsIplQ6a0OUhNORqUfSJd/yq4DF0lNUnn06NEw9swzz4SxF154IYxFI+mWLl0abjNv3rwwlhpxePLkyTAWlRVzHPWmM7tIJpTsIplQsotkQskukgklu0gmdDV+EKKr7qkr7qlYI5aGKjNY54033ghjq1fH63+89tprYWz+/PlV26+77rpwm46OjjCWqnikrqxHv3dqm9TfbDjTmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTAxkRZjpwPeByVRWgFnp7t8xs4uBHwEzqKwK8xl3P9i4rrZeVJJJlWpSJZ5ULDWvWqqMFsUOHoz/NC+++GKpWKr/V1111aDaAUaPHh3GUstQpQbJRMfjfC2vpQzkzN4LfMXd5wIfBr5kZnOBe4A17j4HWFP8LCJDVL/J7u573H1DcfsI8AowDbgFeLC424PArQ3qo4jUwaDes5vZDGARsA6Y3Gcl171UXuaLyBA14GQ3s3HAT4Avu/vhvjGvvAGq+ibIzFaYWZeZdfX09NTUWREpb0DJbmYjqST6D9390aJ5n5lNLeJTgf3VtnX3le7e6e6d7e3t9eiziJTQb7Jb5ZLr/cAr7v6tPqHHgTuK23cAP61/90SkXgYy6u0a4PPAJjPrLtq+CnwdeNjM7gR2Ap9pSA+HgVQJKhqFBumRXO+8804YmzBhQhiLRsS9/PLL4TaPPPJIGHv11VfD2NVXXx3Grr322qrts2bNCrcpOy9cmRFsZUfRDWf9Jru7PwtEv/0N9e2OiDSKPkEnkgklu0gmlOwimVCyi2RCyS6SCU04OQipkWhlpMpyZZc7ikpsqfLamjVrwlhbW1sY+8QnPhHGrrnmmqrt48ePD7dJLUOVmpyzzFJZZUcqDmc6s4tkQskukgklu0gmlOwimVCyi2RCyS6SCZXeBiFVKisjtWZbavLFVAnw9ddfr9qemjjyyJEjYezGG28MYx/72MfC2CWXXFK1PTXaLKXs2nc5TiwZ0ZldJBNKdpFMKNlFMqFkF8mEkl0kE02/Gj+cr45GAy5SAydSV4pTsdRAmP37q07kC8CmTZuqtkdX6QFmzJgRxj71qU+FsdQcdKNGjaranppbL9oG0sf4xIkTYSz6m6Wu7g/n52iKzuwimVCyi2RCyS6SCSW7SCaU7CKZULKLZKLf0puZTQe+T2VJZgdWuvt3zOxe4C7gzNKsX3X31Y3q6PkoVf6J5k4D+M1vfhPGfvnLX1ZtT83vdvvtt4ex5cuXh7EpU6aEsagclhoIU2Yuuf5Ex7hsuXQ4G0idvRf4irtvMLM2YL2ZPVnEvu3u/9q47olIvQxkrbc9wJ7i9hEzewWY1uiOiUh9Deo9u5nNABYB64qmu81so5k9YGbx0qIi0nIDTnYzGwf8BPiyux8G7gNmAQupnPm/GWy3wsy6zKyrp6en2l1EpAkGlOxmNpJKov/Q3R8FcPd97n7K3U8D3wWWVNvW3Ve6e6e7d7a3t9er3yIySP0mu1UuW94PvOLu3+rTPrXP3T4NbK5/90SkXgZyNf4a4PPAJjPrLtq+CtxmZguplON2AF8YyA6H89I6Ufmn7Lxq73//+8PY3r17w9hTTz0Vxrq7u6u2X3nlleE2d911Vxi77LLLwtixY8fCWGTs2LFhLFVeS8VS8/VFz7fU32w4P0dTBnI1/lmg2m+vmrrIMKJP0IlkQskukgklu0gmlOwimVCyi2RCyz8NQlSSSS0LlZo4MjX54hNPPBHGNm+OP9IwbVr1YQuf/OQnw22uuuqqMFZ2dFg0eWTq8VLltdRxTI0ePHnyZNX21KSS9V7ma6jQmV0kE0p2kUwo2UUyoWQXyYSSXSQTSnaRTKj0NghlSm+pstCbb74ZxtauXRvGtmzZEsY6Ojqqtk+dOrVqO6TLUPUuUaVGm6WOVSpWdj29iEpvIjKsKdlFMqFkF8mEkl0kE0p2kUwo2UUy0fTSW6qUM9RF5Z9UOenw4cNhLDV6bcOGDWHs4MGDYWzBggVV21Olt9Ros2j0GqRLVGUm4Uyt9ZYqoaWeU6n+R8pOIDrU6cwukgklu0gmlOwimVCyi2RCyS6SiX6vxpvZGGAtMLq4/4/d/WtmNhNYBUwE1gOfd/fqE36d51JXilNX43fu3BnGUivetrW1hbFomafZs2eH25S9qp7arkzVpRGDXaLlplL9O3HiRBgbzgZyZj8BXO/uV1NZnnm5mX0Y+AbwbXefDRwE7mxYL0WkZv0mu1ccLX4cWXw5cD3w46L9QeDWRnRQROpjoOuzjyhWcN0PPAlsAw65e29xl11A9TmMRWRIGFCyu/spd18IXAosAf5koDswsxVm1mVmXan3oSLSWIO6Gu/uh4CngaXAB8zszAW+S4HdwTYr3b3T3Tvb29tr6auI1KDfZDezdjP7QHF7LPBx4BUqSf/nxd3uAH7aoD6KSB0MZCDMVOBBMxtB5Z/Dw+7+v2b2W2CVmf0z8CJw/0B2mFr+Z6iL+j5mzJhwm8mTJ4exOXPmhLEpU6aEsdSAkUWLFlVtnzlzZrhNyvHjx0v1IyrL9fb2Vm2HdAkttfxTqox29OjRqu1l5w0czvpNdnffCLznGeTu26m8fxeRYeD8/BcmIu+hZBfJhJJdJBNKdpFMKNlFMmHNnBPOzHqAM0O9JgEHmrbzmPpxNvXjbMOtHx3uXvXTa01N9rN2bNbl7p0t2bn6oX5k2A+9jBfJhJJdJBOtTPaVLdx3X+rH2dSPs503/WjZe3YRaS69jBfJhJJdJBMtSXYzW25mvzOzrWZ2Tyv6UPRjh5ltMrNuM+tq4n4fMLP9Zra5T9vFZvakmW0pvk9oUT/uNbPdxTHpNrObm9CP6Wb2tJn91sxeNrO/LdqbekwS/WjqMTGzMWb2vJm9VPTjH4v2mWa2rsibH5nZ4Bayc/emfgEjqMxhdzkwCngJmNvsfhR92QFMasF+PwosBjb3afsX4J7i9j3AN1rUj3uBv2vy8ZgKLC5utwGvAXObfUwS/WjqMQEMGFfcHgmsAz4MPAx8rmj/T+CLg3ncVpzZlwBb3X27V+aZXwXc0oJ+tIy7rwXePKf5Fiqz9EKTZusN+tF07r7H3TcUt49QmQlpGk0+Jol+NJVX1H1G51Yk+zTgD31+buXMtA48YWbrzWxFi/pwxmR331Pc3gvEU9w03t1mtrF4md/wtxN9mdkMKpOlrKOFx+ScfkCTj0kjZnTO/QLdMndfDNwEfMnMPtrqDkHlPzuVf0StcB8wi8qCIHuAbzZrx2Y2DvgJ8GV3P2spnWYekyr9aPox8RpmdI60Itl3A9P7/BzOTNto7r67+L4feIzWTrO1z8ymAhTf97eiE+6+r3iinQa+S5OOiZmNpJJgP3T3R4vmph+Tav1o1TEp9n2IQc7oHGlFsr8AzCmuLI4CPgc83uxOmNmFZtZ25jZwI7A5vVVDPU5lll5o4Wy9Z5Kr8GmacEysMpPn/cAr7v6tPqGmHpOoH80+Jg2b0blZVxjPudp4M5UrnduAv29RHy6nUgl4CXi5mf0AHqLycvBdKu+97qSyQOYaYAvwFHBxi/rxA2ATsJFKsk1tQj+WUXmJvhHoLr5ubvYxSfSjqccEWEBlxuaNVP6x/EOf5+zzwFbgEWD0YB5XH5cVyUTuF+hEsqFkF8mEkl0kE0p2kUwo2UUyoWQXyYSSXSQT/weYte4bQI3fgAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEICAYAAACZA4KlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAASZUlEQVR4nO3df4xdZZ3H8feH2h9AK6XbsU5K6ZS2AmWBgYxdQQS2VoP800o2YmOUjWRrNpBogskSN1nZxGR1s2rchLgpP2I1rsiKRlzrQiEutEJopy5Oi0VoaRFK2xki/QHdpbT97h/3NE7rec7M3J/tPJ9XMpk7z/eee74c+plz5zz3nKOIwMzGvzM63YCZtYfDbpYJh90sEw67WSYcdrNMOOxmmXDYzTLhsNuYSfpzSY9Iel3Sn3xQQ9J/S/o/SW8WX7/rRJ92Iofd6vEO8CBwa8Vzbo+IqcXXhW3qyyo47OOMpJ2SvihpQNJ+ST+UNKWZ64iI30XEfcBzzXxday2HfXz6BHADMA+4DPjrsidJukbSvoqvaxro4Z+Kt/m/knR9A69jTfKuTjdgLfGvEfEagKSfAb1lT4qI9cD0Fqz/74DfAoeBTwI/k9QbEdtbsC4bJe/Zx6c9wx4fAqa2c+UR8UxEHIyItyNiNfAr4MZ29mB/ymHPmKQPDTtiXvb1oSatKgA16bWsTn4bn7GIWEcde31JAiYDk4qfp9ReLt6WNB34C+AJ4AhwM3At8PkmtW11ctitHnOBHcN+/l/gZaAHmAh8BbgIOAo8DyyPiBfa3KOdRL54hVke/De7WSYcdrNMOOxmmXDYzTLR1qPxM2fOjJ6ennau0iwrO3fu5PXXXy/9TENDYZd0A/AtYAJwb0R8ter5PT09bNy4sZFVmlmF97///cla3W/jJU0A7gY+BiwCVkhaVO/rmVlrNfI3+2JgW0S8FBGHgQeAZc1py8yarZGwzwZeGfbzq8XYCSStlNQvqX9oaKiB1ZlZI1p+ND4iVkVEX0T0dXV1tXp1ZpbQSNh3AXOG/XxeMWZmp6BGwr4RWChpnqRJ1C5S8HBz2jKzZqt76i0ijki6HXiE2tTb/RHha5KZnaIammePiDXAmib1YmYt5I/LmmXCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2XCYTfLhMNulgmH3SwTDrtZJhx2s0w47GaZcNjNMuGwm2WioTvCSNoJHASOAkcioq8ZTZlZ8zUU9sJfRsTrTXgdM2shv403y0SjYQ/gUUmbJK0se4KklZL6JfUPDQ01uDozq1ejYb8mIq4EPgbcJunak58QEasioi8i+rq6uhpcnZnVq6GwR8Su4vsg8BNgcTOaMrPmqzvsks6WNO34Y+CjwJZmNWZmzdXI0fhZwE8kHX+df4+I/2pKV2bWdHWHPSJeAi5vYi9m1kKeejPLhMNulgmH3SwTDrtZJprx2fhsHDx4sHT8wIEDyWXOOCP9+7S7u7uuPooZkFJvvPFG6fiuXbvqWlfqvxnglVdeGfNyx44dSy5TtR137NiRrFW95mWXXVY6vmTJkuQyCxcuTNZOZ96zm2XCYTfLhMNulgmH3SwTDrtZJnw0fgy2b99eOr5u3brkMnv27EnWFi1a1HBPJ3vttddKx1O9j2T//v3JWtUR8tQsxNSpU5PL7Nu3L1nbtGlTshYRydqyZctKxy+++OLkMj4ab2anNYfdLBMOu1kmHHazTDjsZplw2M0y4am3MRgcHCwdr5p6e/TRR5O1w4cPN9zTySZPnlw6XjXlNXHixGSt6kSeqtrcuXPHNA4wffr0ZG3btm3J2qFDh5K1888/v3R81qxZyWXGK+/ZzTLhsJtlwmE3y4TDbpYJh90sEw67WSY89TYGS5cuLR2fM2dOcplLL700WVu/fn2yduTIkdE3NszixeW320ud/QUwf/78utZVdbbZWWedVTpeNd34xBNPJGtV2+q9731vsnbdddeVjr/vfe9LLjNejbhnl3S/pEFJW4aNzZC0VtKLxfdzW9ummTVqNG/jvwPccNLYncDjEbEQeLz42cxOYSOGPSKeBP5w0vAyYHXxeDWwvLltmVmz1XuAblZE7C4e76F2R9dSklZK6pfUPzQ0VOfqzKxRDR+Nj9pRmuSRmohYFRF9EdHX1dXV6OrMrE71hn2vpG6A4nv5GSJmdsqod+rtYeAW4KvF9582raNTWOosrwULFiSXue2225K1z372s3X1Uc+U19lnn51cZtKkSXX1USW1rQYGBpLL3Hvvvcla1cUoly9fnqylptiqbqE1Xo1m6u0HwNPAhZJelXQrtZB/RNKLwNLiZzM7hY24Z4+IFYnSh5vci5m1kD8ua5YJh90sEw67WSYcdrNM+Ky3Jqi6YOOMGTPa2MmpI3URyOeffz65TH9/f7L2zjvvJGuXXHJJsvae97yndLxq6q1qavN05j27WSYcdrNMOOxmmXDYzTLhsJtlwmE3y4Sn3k5D9UwbtXuqafv27aXjVffFO3DgQLLW29ubrF1++eXJ2rvf/e5kLTfes5tlwmE3y4TDbpYJh90sEw67WSZ8NP40VM/R81Ycca+6RVXqpJbHHnssuczkyZOTtZtuuilZu+iii5K1M888s3R8vJ7sUsV7drNMOOxmmXDYzTLhsJtlwmE3y4TDbpYJT71Z3fbs2ZOsbdiwoXT8hRdeSC5zwQUXJGtLlixJ1qZNm5as5TjFljKa2z/dL2lQ0pZhY3dJ2iXp2eLrxta2aWaNGs3b+O8AN5SMfzMieouvNc1ty8yabcSwR8STwB/a0IuZtVAjB+hulzRQvM0/N/UkSSsl9UvqHxoaamB1ZtaIesP+bWA+0AvsBr6eemJErIqIvojo6+rqqnN1ZtaousIeEXsj4mhEHAPuARY3ty0za7a6pt4kdUfE7uLHjwNbqp5v49OaNenjsuvXry8d7+7uTi7zqU99KlmrusVT1dly9kcjhl3SD4DrgZmSXgW+DFwvqRcIYCfwuda1aGbNMGLYI2JFyfB9LejFzFrIH5c1y4TDbpYJh90sEw67WSZ81ptVnhlWdUumjRs3Jms7duwoHe/p6UkuU3UbpylTpiRrNjres5tlwmE3y4TDbpYJh90sEw67WSYcdrNMeOptnJFUOl41vVZ1z7Zf/OIXyVrqfm6Qvsfa4sXps6GvuOKKZO2MM7xfapS3oFkmHHazTDjsZplw2M0y4bCbZcJH48eZ1FH3o0ePJpcZHBxM1h566KFkrepWTpdeemnp+NKlS5PLzJ07N1mzxnnPbpYJh90sEw67WSYcdrNMOOxmmXDYzTIxmjvCzAG+C8yidgeYVRHxLUkzgB8CPdTuCvOJiHijda3acamTXSA99fbWW28ll3nqqaeSteeeey5ZqzqBZtGiRaXjvb29yWWstUazZz8C3BERi4APALdJWgTcCTweEQuBx4ufzewUNWLYI2J3RPy6eHwQ2ArMBpYBq4unrQaWt6hHM2uCMf3NLqkHuAJ4Bpg17E6ue6i9zTezU9Sowy5pKvAQ8IWIOOFi4lH7Q7H0j0VJKyX1S+ofGhpqqFkzq9+owi5pIrWgfz8iflwM75XUXdS7gdIPWEfEqojoi4i+rq6uZvRsZnUYMeyqHfq9D9gaEd8YVnoYuKV4fAvw0+a3Z2bNMpqz3j4IfBrYLOnZYuxLwFeBByXdCrwMfKIlHdqfqLqeXGpa7ve//31ymbvvvjtZS93GCWD+/PnJ2nXXXVc6fuGFFyaXsdYaMewRsR5ITex+uLntmFmr+BN0Zplw2M0y4bCbZcJhN8uEw26WCV9wcgyqzjZLqZomq3ddVa85MDBQOn7PPfckl9m0aVOy9vbbbydrK1asSNauv/760vEJEyYkl2n3tsqN9+xmmXDYzTLhsJtlwmE3y4TDbpYJh90sE556G4N2TuPUc2YbpM9ue/rpp5PLHDp0KFm7+uqrk7XU9BrAeeedVzreim3o6bXR8Z7dLBMOu1kmHHazTDjsZplw2M0y4aPxLdaKkzT27t2brG3evLl0/OWXX04ukzpyDvCZz3wmWbv44ouTtXe9q/yflk926Rzv2c0y4bCbZcJhN8uEw26WCYfdLBMOu1kmRpx6kzQH+C61WzIHsCoiviXpLuBvgOO3Zv1SRKxpVaOnq3qnhY4dO5asrVu3Lllbs6b8f0HVteSWL1+erN10003JWtWNOlP/3fVOoXl6rXGjmWc/AtwREb+WNA3YJGltUftmRPxL69ozs2YZzb3edgO7i8cHJW0FZre6MTNrrjH9zS6pB7gCeKYYul3SgKT7JZ3b7ObMrHlGHXZJU4GHgC9ExAHg28B8oJfanv/rieVWSuqX1D80NFT2FDNrg1GFXdJEakH/fkT8GCAi9kbE0Yg4BtwDLC5bNiJWRURfRPRVHdAxs9YaMeyqHT69D9gaEd8YNt497GkfB7Y0vz0za5bRHI3/IPBpYLOkZ4uxLwErJPVSm47bCXyuBf2Na1XTSfv370/WHnnkkWQtdSunSy65JLnMHXfckayde64PxYwXozkavx4omxz1nLrZacSfoDPLhMNulgmH3SwTDrtZJhx2s0z4gpMtVnWW19GjR5O1tWvXJmsbN25M1s4///zS8Ztvvjm5zPz585O11IUjob5bVLX77LVTpY9TgffsZplw2M0y4bCbZcJhN8uEw26WCYfdLBOeemuxqgtHvvnmm8naz3/+82Rt+/btydqCBQtKx6vOXqt3eq1Ksy84Wa8cp9hSvGc3y4TDbpYJh90sEw67WSYcdrNMOOxmmfDU2xhUTRulvPXWW8naU089laxt2LAhWauasjvnnHNKx2fNmpVcpt3TYc12uvffLt6zm2XCYTfLhMNulgmH3SwTDrtZJkY8Gi9pCvAkMLl4/o8i4suS5gEPAH8GbAI+HRGHW9lsp9VzZPfw4fQm2bp1a7K2b9++ZG3q1KnJ2sKFC0vHq27/1Ioj1u289puPuI/OaPbsbwNLIuJyardnvkHSB4CvAd+MiAXAG8CtLevSzBo2Ytij5vjE7sTiK4AlwI+K8dXA8lY0aGbNMdr7s08o7uA6CKwFtgP7IuJI8ZRXgdkt6dDMmmJUYY+IoxHRC5wHLAYuGu0KJK2U1C+pf2hoqL4uzaxhYzoaHxH7gF8CVwHTJR0/wHcesCuxzKqI6IuIvq6urkZ6NbMGjBh2SV2SphePzwQ+AmylFvq/Kp52C/DTFvVoZk0wmhNhuoHVkiZQ++XwYET8p6TfAg9I+grwP8B9LezztHXWWWcla1dffXWyNnPmzGRt9uz04ZGrrrqqdHzevHnJZeqdumr2CSitOKHFt3/6oxHDHhEDwBUl4y9R+/vdzE4D/gSdWSYcdrNMOOxmmXDYzTLhsJtlQu2cgpA0BLxc/DgTeL1tK09zHydyHyc63fqYGxGln15ra9hPWLHUHxF9HVm5+3AfGfbht/FmmXDYzTLRybCv6uC6h3MfJ3IfJxo3fXTsb3Yzay+/jTfLhMNulomOhF3SDZJ+J2mbpDs70UPRx05JmyU9K6m/jeu9X9KgpC3DxmZIWivpxeL7uR3q4y5Ju4pt8qykG9vQxxxJv5T0W0nPSfp8Md7WbVLRR1u3iaQpkjZI+k3Rxz8W4/MkPVPk5oeSJo3phSOirV/ABGrXsLsAmAT8BljU7j6KXnYCMzuw3muBK4Etw8b+GbizeHwn8LUO9XEX8MU2b49u4Mri8TTgBWBRu7dJRR9t3SaAgKnF44nAM8AHgAeBTxbj/wb87VhetxN79sXAtoh4KWrXmX8AWNaBPjomIp4E/nDS8DJqV+mFNl2tN9FH20XE7oj4dfH4ILUrIc2mzdukoo+2ipqmX9G5E2GfDbwy7OdOXpk2gEclbZK0skM9HDcrInYXj/cA6Xsst97tkgaKt/kt/3NiOEk91C6W8gwd3CYn9QFt3iatuKJz7gforomIK4GPAbdJurbTDUHtNzu1X0Sd8G1gPrUbguwGvt6uFUuaCjwEfCEiDgyvtXOblPTR9m0SDVzROaUTYd8FzBn2c/LKtK0WEbuK74PAT+jsZbb2SuoGKL4PdqKJiNhb/EM7BtxDm7aJpInUAvb9iPhxMdz2bVLWR6e2SbHufYzxis4pnQj7RmBhcWRxEvBJ4OF2NyHpbEnTjj8GPgpsqV6qpR6mdpVe6ODVeo+Hq/Bx2rBNVLsq5H3A1oj4xrBSW7dJqo92b5OWXdG5XUcYTzraeCO1I53bgb/vUA8XUJsJ+A3wXDv7AH5A7e3gO9T+9rqV2g0yHwdeBB4DZnSoj+8Bm4EBamHrbkMf11B7iz4APFt83djubVLRR1u3CXAZtSs2D1D7xfIPw/7NbgC2Af8BTB7L6/rjsmaZyP0AnVk2HHazTDjsZplw2M0y4bCbZcJhN8uEw26Wif8HhJAtV2Fb+ucAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 然后我们看看经典奇异值的分解效果\n",
    "U, sigma, V = np.linalg.svd(imgmat)\n",
    "\n",
    "for i in range(5, 16, 5):\n",
    "    reconstimg = np.matrix(U[:, :i]) * np.diag(sigma[:i]) * np.matrix(V[:i, :])\n",
    "    plt.imshow(reconstimg, cmap='gray')\n",
    "    title = \"n = %s\" % i\n",
    "    plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:18.966739Z",
     "start_time": "2021-03-09T03:50:18.950606Z"
    }
   },
   "outputs": [],
   "source": [
    "# 然后我们再来看看量子版本的分解效果：\n",
    "\n",
    "# 超参数设置\n",
    "N = 5           # 量子比特数量\n",
    "T = 8           # 设置想要学习的阶数\n",
    "ITR = 200       # 迭代次数\n",
    "LR = 0.02        # 学习速率\n",
    "SEED = 14       # 随机数种子\n",
    "\n",
    "# 设置等差的学习权重\n",
    "weight = np.arange(2 * T, 0, -2).astype('complex128')\n",
    "\n",
    "\n",
    "def Mat_generator():\n",
    "    imgmat = np.array(list(img.getdata(band=0)), float)\n",
    "    imgmat.shape = (img.size[1], img.size[0])\n",
    "    lenna = np.matrix(imgmat)\n",
    "    return lenna.astype('complex128')\n",
    "\n",
    "M_err = Mat_generator()\n",
    "U, D, V_dagger = np.linalg.svd(Mat_generator(), full_matrices=True)\n",
    "\n",
    "# 设置电路参数\n",
    "cir_depth = 40                           # 电路深度\n",
    "block_len = 1                            # 每个模组的长度"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:50:19.553293Z",
     "start_time": "2021-03-09T03:50:19.531180Z"
    }
   },
   "outputs": [],
   "source": [
    "# 重新定义量子神经网络\n",
    "def U_theta():\n",
    "\n",
    "    # 用 Circuit 初始化网络\n",
    "    cir = Circuit(N)\n",
    "    \n",
    "    # 搭建层级结构：\n",
    "    for _ in range(cir_depth):\n",
    "        cir.ry()\n",
    "        cir.cnot()\n",
    "\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-09T03:59:58.649381Z",
     "start_time": "2021-03-09T03:54:30.126561Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 0 loss: -4291.2284\n",
      "iter: 10 loss: -103444.3877\n",
      "iter: 20 loss: -126821.2327\n",
      "iter: 30 loss: -135986.5470\n",
      "iter: 40 loss: -141457.8124\n",
      "iter: 50 loss: -145045.8896\n",
      "iter: 60 loss: -147249.5471\n",
      "iter: 70 loss: -149020.3769\n",
      "iter: 80 loss: -150570.5096\n",
      "iter: 90 loss: -151898.1960\n",
      "iter: 100 loss: -152989.3110\n",
      "iter: 110 loss: -153869.6654\n",
      "iter: 120 loss: -154547.8973\n",
      "iter: 130 loss: -155069.0768\n",
      "iter: 140 loss: -155459.7348\n",
      "iter: 150 loss: -155744.3048\n",
      "iter: 160 loss: -155957.9127\n",
      "iter: 170 loss: -156124.0155\n",
      "iter: 180 loss: -156257.4063\n",
      "iter: 190 loss: -156367.2717\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1d251eaff70>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD5CAYAAADhukOtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX4ElEQVR4nO2da4yUZZbH/wdEbt2A3VyFVi6LbMiwILRGM2rcMTNxzSRqsjH6wfjBDJPNmKzJ7AfjJqub7Adns2r8sHGDqxln43rZcYxmorvjmknIaHRsEBBEF+Si3TQ09+4GvABnP9RLpiF1/l08VfUW8vx/CaH6OfXUe+qt999d9fzrnMfcHUKIi58xrU5ACFEOErsQmSCxC5EJErsQmSCxC5EJErsQmXBJPZPN7FYATwEYC+Df3f0xdv+Ojg7v6uqqGjt9+nQ9qTQMMwtjY8ZU/93I5jBrkz3nb7/99rzzAICxY8ee1zgAnDp1KoylWrPROWG5p+bBnluUBzv37FjstWY02uKO8ujt7cWhQ4eqBpPFbmZjAfwrgB8C6AXwoZm94e6fRHO6urrw1ltvVY19/fXX4bGii4C9yKkn99JLLw1j48ePP69xgIt2eHg4jA0MDIQxlmNnZ2fV8ba2tqQ8vvnmmzDGxDlx4sSq4yx3lsfJkyfDGHtul1xS/RJn1xt7ztHjAfx6ZNdBBDu/EyZMqDp+2223hXPqeRt/LYDt7r7D3b8B8BKA2+t4PCFEE6lH7HMBfDni595iTAhxAdL0BTozW21mPWbWc/DgwWYfTggRUI/Y+wCMXG2bV4ydhbuvcfdud++OPk8KIZpPPWL/EMBiM1tgZpcCuBvAG41JSwjRaJJX4939pJk9AOB/ULHennP3LaPNS7GvolXOVDuGwayhCJZ7yvMCgHHjxoUxtrI7NDRUdZyt7KbYjaPFImvrxIkT4Ry2Cp56HqPVc3Y+WCz1umKr+NHryeYk5VDPZHd/E8CbDcpFCNFE9A06ITJBYhciEyR2ITJBYhciEyR2ITKhsWv7NRBZMqwKKYqlWD9AumUXzSu7Yo9ZMlEBSjQOcMuLnStmUaVYrKnnkc2Lcky13tg1x55bSkUfK/5JseX0l12ITJDYhcgEiV2ITJDYhcgEiV2ITCh1Nd7dk1bWU4pTGKkrwilFEGylmxWFsGKXKVOmhLH29vaq4+w5Hz9+PCkPVqyTkgdbfT569GgYGxwcDGOTJ0+uOp7a44+tqrPzwYjamqXkSIuazi8tIcR3FYldiEyQ2IXIBIldiEyQ2IXIBIldiEwovRAmgtkMKbt6sF1aGMxaiWCWHLNxol09AGDatGlhjO2qEuXfjF1OIlsLiC0gZjemWoApvfyYzdeMLbtSt9+KiPJnuesvuxCZILELkQkSuxCZILELkQkSuxCZILELkQl1WW9mtgvAEIBTAE66e/docyKbIaWnVltbWxhL6WkHpNlozO5g9hrLg1V5pVSipW4nxaqoWLVZ9Hqyc8UsRdZDj8Wi/FMrH9l1mlotFz1mSh7seTXCZ/9Ldz/QgMcRQjQRvY0XIhPqFbsD+J2ZrTOz1Y1ISAjRHOp9G3+Du/eZ2UwAb5vZp+6+duQdil8CqwFg7ty5dR5OCJFKXX/Z3b2v+H8AwGsArq1ynzXu3u3u3R0dHfUcTghRB8liN7PJZtZ+5jaAHwHY3KjEhBCNpZ638bMAvFYs9V8C4D/d/b/ZBHcPLQhmGUR2GK3wSWxSmVKBlGrzpVZesWo/VlUWwezGlCabQFxlx2y+1O2fUqzDVPs15fpgeQBp1lsUY9dNstjdfQeA5anzhRDlIutNiEyQ2IXIBIldiEyQ2IXIBIldiEwoteGkmYW2RorFwyrKWPPC1H3lhoeHq47v3bs3nMOq1yZNmhTGWCXXZZddFsaiyrGhoaFwzs6dO8MYy//w4cNhLDqPR44cCef09vaGsa+++iqMzZ8/P4wtX17dMOrq6grnsEaarDKPkVLBxizWSC/a600IIbELkQsSuxCZILELkQkSuxCZUOpqvLuHq5IpWxqxogS2kslWullRRV9fX9Xx9evXh3N27doVxlh9PysHZs87WtFmK+f79+8PY2weyyNa0War+9u2bQtje/bsCWPLli0LY9EqPju/zOVh24qxLbbYKvmxY8eqjrNCqSjGVv31l12ITJDYhcgEiV2ITJDYhcgEiV2ITJDYhciEUq03ILZrmPUW2WGpW/iwYhdmd0RFHJElBwDbt28PY8xqOnAg3mSH2WFRjFlNrAgpKv4BgNmzZ4exmTNnVh3v7OwM57D+aaxIhtlyUZESs2YZ7PpgW0OxedH1yK7TyOajfRnDiBDiokJiFyITJHYhMkFiFyITJHYhMkFiFyITRrXezOw5AD8GMODu3yvGOgC8DGA+gF0A7nL32A8qGDNmDNra2qrGUip8Urd/YrYcq4iLep2xY02fPj2MMcuIVd8tXrw4jEUVYLfccks4h1V5sXM1derUMBb1yWNVgC+//HIY6+/vD2NLliwJY6tWrao63t7eHs5h1xWD2cfs+o6sT1ZFF1UV0v6KYeRP/BLAreeMPQTgHXdfDOCd4mchxAXMqGIv9ls/dM7w7QCeL24/D+COxqYlhGg0qZ/ZZ7n7mfdVe1HZ0VUIcQFT9wKdVz7ghB9yzGy1mfWYWc/BgwfrPZwQIpFUse8zszkAUPw/EN3R3de4e7e7d7PvRQshmkuq2N8AcF9x+z4ArzcmHSFEs6jFensRwM0ApptZL4BHADwG4BUzux/AbgB31XIwM8O4ceOqxth2TZHNcOjQueuGfyKy+ABukZw4cSKMRTZa9JwAbl1df/31YYy9C2LVVdF5ZNbPlClTwhhrsMi2SYqssvfffz+c89lnnyUda9GiRWEsskvZc2b2K7PlmFXGjhdV4LFrOKr4ZNbbqGJ393uCUGzcCiEuOPQNOiEyQWIXIhMkdiEyQWIXIhMkdiEyofSGkxGs6WFUATZp0qRwDtuHjMWY1RTZaOzxmM3HLDT2mMyWixo9Mnsw2g8N4HbSl19+GcYiW5RZm+xY7Dyy55ZSwcZsypTrA+BVjClzmD0Yob/sQmSCxC5EJkjsQmSCxC5EJkjsQmSCxC5EJpRqvZlZaDexvbeiOYODg+GcqOEhkL7XWwSryGJ5RDYZwKukjh49GsaiajNmGS1cuDCMMeuK2Wiffvpp1fF169aFc9gedmxfuQULFoSxyKZMqSgbLcasMnaNRM+bWctJTVjDiBDiokJiFyITJHYhMkFiFyITJHYhMqHU1fhTp05haGioaiylEIYVA7CiBFZkkgJbwZ82bVoYY8Ud0XkCuAsRrbrPmDEjKQ+24r5jx44wtm3btqrjR44cCed0dHSEsauvvjqMzZs3L4xFTg57zdiKe+q2YseOHQtj0bXPcmTHitBfdiEyQWIXIhMkdiEyQWIXIhMkdiEyQWIXIhNq2f7pOQA/BjDg7t8rxh4F8BMA+4u7Pezub9ZywOiL+swOmzBhQtVxZk8xW44di/Vji/JItQDb29vDGLPDWNFQZPHs2bMnnDN16tQwxs4V66EXWWxsJ9/LL788jEVbbwG8J19UTEK3SSLPi1mR7HVJuR5TrNl6C2F+CeDWKuNPuvuK4l9NQhdCtI5Rxe7uawHEOygKIb4T1POZ/QEz22Rmz5lZXLQthLggSBX70wAWAVgBoB/A49EdzWy1mfWYWQ/bYlkI0VySxO7u+9z9lLufBvAMgGvJfde4e7e7d7PvPgshmkuS2M1szogf7wSwuTHpCCGaRS3W24sAbgYw3cx6ATwC4GYzWwHAAewC8NNaDxjZDCl2R0ql3GgxVvEUVSExe4rZa2zbImYrMvsnOo8sj+Hh4TC2b9++MLZx48YwtnPnzjAWccUVV4SxVatWhbErr7wyjKX0FGTbUKXCXrPo+knplcist1HF7u73VBl+drR5QogLC32DTohMkNiFyASJXYhMkNiFyASJXYhMKLXhJBDbZazCJ6pEY3OYBcFgFmBkkTB7h9lyzI5h3zY8fvx4GIu2NWI2H9tOau/evWGsr68vjEWvGbMA2Zeu2HZY7DVjNmsK7FjsmmOxyO6lFWyBLUebZYYRIcRFhcQuRCZI7EJkgsQuRCZI7EJkgsQuRCaUar25e1hRlGKfTJ48OZzDmv+lVEIBsRXCHo/ZJ6n7jc2cOTOMRfYVy2P37t1hLNqzbbR5Uf6zZ88O58ydOzeMXXZZ3AyJVYdF1w6zRNnjMVhzUWYBRo1MWfVddKx6G04KIS4CJHYhMkFiFyITJHYhMkFiFyITSl2NHzNmTLgNTkqhACskSXm80eZFq+fs8QYHB8MYK+SJVmiBuNgFiFeZv/jii3DOpk2bwtiGDRvCGCvIiVbdlyxZEs5hq/FsZTrVDYlgrycjdUup6DVLvU4j9JddiEyQ2IXIBIldiEyQ2IXIBIldiEyQ2IXIhFq2f+oC8CsAs1DZ7mmNuz9lZh0AXgYwH5UtoO5y98PssVghDLOhoiKCzs7OcE7qjrGsV1tE1G8N4BYas65YIU9kX7IYy2NgYCCMHThwIIyx3nULFy6sOr506dJwDtvGaerUqWGM2ZvReWRbhzGbjF2nzHpjr3VU0MWs5SiPenvQnQTwc3dfCuA6AD8zs6UAHgLwjrsvBvBO8bMQ4gJlVLG7e7+7ry9uDwHYCmAugNsBPF/c7XkAdzQpRyFEAzivz+xmNh/A1QA+ADDL3fuL0F5U3uYLIS5Qaha7mbUBeBXAg+5+1ockr3x3r+r398xstZn1mFnPwYMH60pWCJFOTWI3s3GoCP0Fd/9NMbzPzOYU8TkAqq7yuPsad+929262oCaEaC6jit0qy3vPAtjq7k+MCL0B4L7i9n0AXm98ekKIRlFL1dv3AdwL4GMz21CMPQzgMQCvmNn9AHYDuKumAwa2RorlxbYSYtYbs12YtRLZhs3od5fany7KkdlrH330URjbsmVLGGOvWVTBtnz58nDO9OnTwxizrpgFGNm2zL5kveTGjx8fxth1dezYsfOex/JIsd5GFbu7/wFA9Ai3jDZfCHFhoG/QCZEJErsQmSCxC5EJErsQmSCxC5EJpTacBOJGecPDw+GcqGJr//794RxWgcQsI5ZHlPukSZPCOUNDQ2GMbWk0Z86cMMaqod59992q4y+++GI4Z8eOHWGMncebbropjEUWG2uWyThy5EgYS6k6pBYVqXpjNiuzythWZZFdymy+aI62fxJCSOxC5ILELkQmSOxCZILELkQmSOxCZMIFY70xyyCq8GEWFLPXmLXCbJyoOin1WKxKKtr/CwB6e3vD2LZt26qO9/f3Vx0HgI6OjjDGrDK2b1tkebEKNVbpx2IpDSLZHPZ6steFxZg9Gz03lkdUzSfrTQghsQuRCxK7EJkgsQuRCRK7EJlQ6mq8uyf1ZIuKMdiWRgy2Cp5SQMOeE+tpx2AFNHv37g1jUcEIe84zZ84MY/PmzQtjCxYsCGMzZsyoOs5WmNl5ZEUmLNZoWAENizEiZ4C9ZtGqe73bPwkhLgIkdiEyQWIXIhMkdiEyQWIXIhMkdiEyYVTrzcy6APwKlS2ZHcAad3/KzB4F8BMAZxrBPezub7LHYtYbKyKI7CtWtMIKD9ixmI0TzWNFFWyLKtZjbHBwMIy99957YWzt2rVVx9n2TzfeeGMYu+aaa8LYypUrw1i0XRPrG8isVFaQM23atDAWWansdWZ95pg1y+ax/CO7lD2vlB50tfjsJwH83N3Xm1k7gHVm9nYRe9Ld/6WGxxBCtJha9nrrB9Bf3B4ys60Aqu/aJ4S4YDmvz+xmNh/A1QA+KIYeMLNNZvacmcV9kYUQLadmsZtZG4BXATzo7oMAngawCMAKVP7yPx7MW21mPWbWc/jw4fozFkIkUZPYzWwcKkJ/wd1/AwDuvs/dT7n7aQDPALi22lx3X+Pu3e7ezTZFEEI0l1HFbpVv1j8LYKu7PzFifOSWJXcC2Nz49IQQjaKW1fjvA7gXwMdmtqEYexjAPWa2AhU7bheAn9aTCKvwiWLUZiB2WGqVVJRH6rGiPmIAcPDgwTB24MCBMBZZMp2dneGcRYsWhbGrrroqjM2ePTuM7d69u+r4sWPHwjmsYotZqcwOY7GI1NcztSIuspZZhSCz+SJqWY3/A4BqmVJPXQhxYaFv0AmRCRK7EJkgsQuRCRK7EJkgsQuRCaU2nDSz0E5IsUhSG/ylbDXF5qVaecxC27hxYxiLbC0AmDhxYtXxFStWhHOWLVsWxpi9xhpERlV7zG5ksOsjpTEjm8OeF8uDXTvsGkm59pPmnPcMIcR3EoldiEyQ2IXIBIldiEyQ2IXIBIldiEwo3XqLLA9W1XTixImq46w+nlk8rOEkyyOq2IryA3jlUm9vbxiLGkcCwOeffx7Goio11jiSNTZkNuWePXvCWGQ1zZo1K5wzefLkMJZaARa91pFFyeaMNo81zGTNUaMqO3YNR+dDe70JISR2IXJBYhciEyR2ITJBYhciEyR2ITKhVOvN3ZMa5UVVSNEeWQCvamLWStSwEYjtJNagkFlXKZYRwKurIvsqZQ7A96Njll2UP6vWmjJlShhjVWPsHEfXDrNLWdXb0NBQGGPXQUq1X8rzYnP0l12ITJDYhcgEiV2ITJDYhcgEiV2ITBh1Nd7MJgBYC2B8cf9fu/sjZrYAwEsAOgGsA3Cvu8dL2QXRKjlb9Y1WMlnRCis8YLBV62j1nK3sMtgWT6zPHFtx7erqqjrOetCx1XjmTrD8o1V3VrzEVuqZ88K2lIpW8dn1lrqtGCuSYddIdP7ZuW9WIczXAH7g7stR2Z75VjO7DsAvADzp7n8G4DCA+2t4LCFEixhV7F5huPhxXPHPAfwAwK+L8ecB3NGMBIUQjaHW/dnHFju4DgB4G8DnAI64+xlnvxfA3KZkKIRoCDWJ3d1PufsKAPMAXAvgz2s9gJmtNrMeM+s5dOhQWpZCiLo5r9V4dz8C4PcArgcwzczOrFbMA9AXzFnj7t3u3t3R0VFPrkKIOhhV7GY2w8ymFbcnAvghgK2oiP6vi7vdB+D1JuUohGgAtRTCzAHwvJmNReWXwyvu/lsz+wTAS2b2TwA+AvBsLQeMrAFW+BHZJ8wiYaTOi2xD1nuM2UmsH9uSJUvCWLS1EntM1sONFQ0xqyml8IPZZAxWCMNstJSCHHYtMmsrdV5KUUvK1mejit3dNwG4usr4DlQ+vwshvgPoG3RCZILELkQmSOxCZILELkQmSOxCZIKl2lBJBzPbD+BMOdd0AAdKO3iM8jgb5XE237U8rnT3GdUCpYr9rAOb9bh7d0sOrjyUR4Z56G28EJkgsQuRCa0U+5oWHnskyuNslMfZXDR5tOwzuxCiXPQ2XohMaInYzexWM/vMzLab2UOtyKHIY5eZfWxmG8ysp8TjPmdmA2a2ecRYh5m9bWbbiv/jzozNzeNRM+srzskGM7uthDy6zOz3ZvaJmW0xs78txks9JySPUs+JmU0wsz+a2cYij38sxheY2QeFbl42s7g7ajXcvdR/AMai0tZqIYBLAWwEsLTsPIpcdgGY3oLj3gRgJYDNI8b+GcBDxe2HAPyiRXk8CuDvSj4fcwCsLG63A/g/AEvLPickj1LPCQAD0FbcHgfgAwDXAXgFwN3F+L8B+JvzedxW/GW/FsB2d9/hldbTLwG4vQV5tAx3Xwvg3B5dt6PSuBMoqYFnkEfpuHu/u68vbg+h0hxlLko+JySPUvEKDW/y2gqxzwXw5YifW9ms0gH8zszWmdnqFuVwhlnu3l/c3gsg7mzRfB4ws03F2/ymf5wYiZnNR6V/wgdo4Tk5Jw+g5HPSjCavuS/Q3eDuKwH8FYCfmdlNrU4IqPxmR+UXUSt4GsAiVPYI6AfweFkHNrM2AK8CeNDdz2rHU+Y5qZJH6efE62jyGtEKsfcBGLltSdisstm4e1/x/wCA19Dazjv7zGwOABT/D7QiCXffV1xopwE8g5LOiZmNQ0VgL7j7b4rh0s9JtTxadU6KYx/BeTZ5jWiF2D8EsLhYWbwUwN0A3ig7CTObbGbtZ24D+BGAzXxWU3kDlcadQAsbeJ4RV8GdKOGcWKWh2rMAtrr7EyNCpZ6TKI+yz0nTmryWtcJ4zmrjbaisdH4O4O9blMNCVJyAjQC2lJkHgBdReTv4LSqfve5HZc+8dwBsA/C/ADpalMd/APgYwCZUxDanhDxuQOUt+iYAG4p/t5V9TkgepZ4TAH+BShPXTaj8YvmHEdfsHwFsB/BfAMafz+PqG3RCZELuC3RCZIPELkQmSOxCZILELkQmSOxCZILELkQmSOxCZILELkQm/D965hTK7FIF7QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 记录优化中间过程\n",
    "loss_list, singular_value_list = [], []\n",
    "U_learned, V_dagger_learned = [], []\n",
    "    \n",
    "net = NET(Mat_generator(), weight)\n",
    "\n",
    "# 一般来说，我们利用 Adam 优化器来获得相对好的收敛\n",
    "# 当然你可以改成 SGD 或者是 RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# 优化循环\n",
    "for itr in range(ITR):\n",
    "\n",
    "    #  前向传播计算损失函数\n",
    "    U, V_dagger, loss, singular_values = net()\n",
    "\n",
    "    # 反向传播极小化损失函数\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # 记录优化中间结果\n",
    "    loss_list.append(loss[0][0].numpy())\n",
    "    singular_value_list.append(singular_values)\n",
    "    \n",
    "    if itr% 10 == 0:\n",
    "        print('iter:', itr,'loss:','%.4f'% loss.numpy()[0])\n",
    "\n",
    "# 记录最后学出的两个酉矩阵    \n",
    "U_learned = U.numpy()\n",
    "V_dagger_learned = V_dagger.numpy()\n",
    "\n",
    "singular_value = singular_value_list[-1]\n",
    "mat = np.matrix(U_learned.real[:, :T]) * np.diag(singular_value[:T])* np.matrix(V_dagger_learned.real[:T, :])\n",
    "\n",
    "reconstimg = mat\n",
    "plt.imshow(reconstimg, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] Wang, X., Song, Z., & Wang, Y. Variational Quantum Singular Value Decomposition. [Quantum, 5, 483 (2021).](https://quantum-journal.org/papers/q-2021-06-29-483/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
